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Introduction 


This unit is devoted to the study of numerical methods for finding approximate 
solutions of differential equations. Only first-order differential equations of the 
form 


y’ = m(x, y), 


with y(x9) given, will be considered but the methods in this unit can be adapted to 
find approximate solutions of higher-order differential equations. (A differential 
equation of any order may be treated as a system of first-order differential 
equations, as described in Section 5 of Unit 6.) 


The primary objective of this unit is to enable you to use numerical methods 
intelligently to obtain approximate solutions. To do this it is necessary to study 
several different methods 

(i) to see how they work, 

(ii) to determine whether one method is better than another, 

(iii) to see when things go wrong. 


All the methods in this text lead to recurrence relations which can be used step by 
step to obtain an approximate solution at equally spaced points Xo, x1, X2,--. . 
As in Unit 18 we denote the values of this approximate solution by Yo, ¥;, ¥)... . 
It is important to compare this solution with the true solution y = y(x). In 
particular we are interested in the difference between the approximate value, Y,,, 
and the true value, y,, at X,. 


Study guide 


This unit is intended to be read sequentially. In Sections 1 and 2 we derive some 
numerical methods and use them to solve simple problems. The television section 
(Section 2) illustrates some of these methods and how they can be used to solve 
the logistic equation. Before viewing the programme it is advisable to have read 
Subsection 1.1 on the trapezoidal method and also the pre-broadcast notes. 


In Section 3 and in the tape session in Section 4 we look at the methods of the 
first two sections to discover why they work, when they go wrong and what sort of 
accuracy you could expect. 


In the remainder of Section 4 we look at Simpson’s method which is more 
accurate than most of the one-step methods but can be difficult to use. 


Section 5 describes a computer package which you can use to solve differential You can use the computer 
equations numerically. One of the TMA questions for this unit may well require package at Summer School. 
you to use this package. 


1 Numerical methods 


1.1. Integration methods 


In Unit 18 you met three integration methods for evaluating the integral 


l= f (x) dx. 
They were Euler’s method, the trapezoidal method and Simpson’s method. In this 
section we will see how the first two of these methods can be used to solve first- 
order differential equations of the form 


d 
y’ = m(x, y). (1) — y’ is short for 


In Section 4 we will examine Simpson’s method. 


If we integrate both sides of Equation (1) between x, and x,,, we have 


{ y dx = | m(x, y) dx. (2) 


r xr 
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The left-hand side of this equation can be integrated to give 


: Y dx =Yyi1—JYy. y, 1s short for y(x,). 


Xr 


Hence we may write Equation (2) as 


Veua Je -{ m(x, y) dx. (3) 


Xr 


Thus, if we can determine the integral on the right-hand side and if we know y,, 
we can compute y,,;, using Equation (3). The difficulty with this approach is that 
we do not know the function y, so it is unlikely that we will be able to determine 
the integral. 


For example, the differential equation 
y= xy with y(0) = 0 
has m(x, y) = xy. 


Equation (3) becomes 


Xr+i1 


Yror =e | xy dx. 


Xr 


Since we do not know y(x) we cannot proceed to determine an exact recurrence 
relation, but we can determine an approximate one if we approximate the integral. 


In the last unit we saw that the simplest method for approximating integrals was 


Euler’s method (see Figure 1) giving 7 
J (x) dx = hf (x,) 
mmcre h= x, — x. f(x,) 


In the same way we can approximate the integral in Equation (3) using Euler’s 
method as 


Xp +1 
Xr+1° h 
[  mex,y) dx ~ hmx,,y,) | | 
- Figure 1. The shaded area is 
: : Euler’s approximation to the 
With this approximation Equation (3) becomes aeea Sa the curve 


: ae = Zz. 5 hm(x,, ¥,). 


(I have used capital letters, Y, and Y,,,, here because Y, will only be an 
approximation to the true solution.) But this recurrence relation is just Euler’s 
method for solving a first-order differential equation, as described in Units 2 and 
18. This new derivation illustrates how integration methods can be used to solve 
differential equations. 


To simplify the notation, we define 
r= m(x,, Y,) 
and write the recurrence relation for Euler’s method as 


Y,+1 — Je r ay. 
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Euler’s method 
1. You are given a differential equation 
y' = m(x, y) 
and a value for y(xo). 


2. To apply Euler’s method, choose a step-size, h, and calculate 
¥y, Yo... rom 


om ; oes Y, + hY, 
where tT; —_ m(X,, ¥)s 


Yo a y(Xo), 
and X,=Xo trh. 


3. Y, is an approximation to y,. 


Another integration method is the trapezoidal method. From the results in Unit 18 
it is clear that this method for approximating the integral by a trapezium (see 
Figure 2) as 


Xr+i1 h 
[ ferde = 511 Ce) + fora] 
is generally more accurate than Euler’s method. We would thus expect to get more 
accurate results if we use this method to solve differential equations. 
If we apply this method to approximate the integral in Equation (3) we have 


Xr+1 h 
[  mtx.y)dx = 5 frm(x p90) + mle 1dr) 


Xr 


This gives rise to the recurrence relation 
h 
Y+1 i 1; + Dy [m(x,, 7) + MUX e441, ) eee 
or, more simply, 
: h / / 
ee ae Y¥,+>5L¥, + Tree (4) 


However there is a difficulty in using the trapezoidal method because Equation (4) 
contains a term in Y},, (= m(x,+1, Y,+1)) on the right-hand side. We do not 
know this value until we have computed Y,,,. Equation (4) would have to be 
solved for Y,,,. For linear equations this can be done fairly easily, as Example 1 
shows. 


Example 1 


Use the trapezoidal method with h = 0.1 to determine the approximate solution 
for x, =rh, Osra hk, of 


y=xt+y with y(0) = 0. 

Compare the answers with the true solution, y = e* — x — 1. 

Solution 

We want to apply the recurrence relation (4). The differential equation tells us that 
m(x,y)=xt+y 

and so Y¥, =x,+ 7, 

and eet 2 Meat Pt te: 


We substitute for Y}, and Y}., in Equation (4) to get 


h 
Fett — i - 73 Lr - Y,) - (Xr41 - ae i)]. 


Figure 2. The shaded area is 
the trapezoidal approximation 
to the area under the curve. 
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To solve for Y,,, we collect up the terms in Y,,, and Y, as 


Yres(1 - 4 = (1 : 4 + =x + X;41). 
ee ge ee) 
1 —h/2 
We can use this recurrence relation to generate the approximate solution. We have 
i) h=0.1, 


Gi) Xx. = ay-+ 7h = 0.1r. 


We can simplify the recurrence relation to 


1.05Y, + 0.05(0.1r + 0.1(r + 1)) 
0.95 


Y,+1 _ 


210Y, + 2r+1 


i.e. | — 190 


The computation is now straightforward, starting with Y) = 0. The following table 
shows how the approximate solution, generated using the recurrence relation for 
the trapezoidal method, compares with the true solution. The results for Euler’s 
method and the Taylor series method of order 2 are also given for comparison. 
(These were obtained in Unit 18.) 


Solutions to y’ = x + y with y(0) = 0.h = 0.1. 


trapezoidal true 
a method solution | Euler’s method 
0.0 0.0000 


Taylor series method 
of order 2 


0 : 
1 | 0.1 0.0000 
2402 0.0100 0.0210 
3 0.0310 0.0492 
~ 0.0641 0.0909 
5 0.1105 0.1474 
6 0.1716 0.2204 
7 0.2487 0.3116 
8 0.3436 0.4228 
9 0.4579 0.5562 
10 0.5937 0.7141 


By comparing columns 3 and 4 we see that the results for the trapezoidal method 
are reasonably accurate although, as expected, the error increases as x increases. 
The approximate value of y at x = 1 is correct to two decimal places. Now 
compare the true solution with the results obtained using Euler’s method. For this 
problem Euler’s method gives much worse results and a much smaller step-size 
would be needed to obtain reasonable accuracy. On the other hand the results for 
the Taylor series method of order 2 are about as accurate as the results for the 
trapezoidal method. 


The manipulation required to use the trapezoidal method to solve the general 
linear differential equation of the form 


y’ = I(x)y + k(x), (5) 


using a step-size h, can be deduced in the same way as in Example 1. The 
recurrence relation (4) for the trapezoidal method is 


h 
144 = Y, + ree + ae 1). (6) 


Using Equation (5) we have 
Y,=1,Y, +k, 
and Yeet 1512 i1 + he. 


This can be expressed as 


ee ok Eee 

ag 85. 41 
which is a linear recurrence 
relation. 


!, is short for I(x,), k, is short 
for k(x,) and so on. 
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Substituting for Y; and Y;,, in Equation (6) gives 
h 
Y+i=¥%t+ alley, +) +(e ita Fhe 0] 


which can be solved for Y,,, as 


Ut ua IY, + (/2Vke + kes) 
aS 1 — (2a 7 


The trapezoidal method for linear differential equations 


1. You are given a linear differential equation 


y’ =I(x)y + k(x) 
and a value for y(xo). 


2. To apply the trapezoidal method, choose a step size h and 
calculate Y,, Y,,... from 


LL (uM 1Y, + (H/2Wke + kro) 
eee a es 


where Yo = y(Xo), and x, = Xo + rh. 


3. Y, is an approximation to y,. 


We have seen two methods in this section: Euler’s method, in which the recurrence 
relation can be applied directly, and the trapezoidal method, which is difficult to 
use except for linear equations. The difference between the two methods, which 
makes the trapezoidal method harder to use, can be described by saying that the 
trapezoidal method is implicit whereas Euler’s is explicit. 


In general any numerical method for solving differential equations in which the 
recurrence relation involves Y’,, or higher derivatives at x,+, is called implicit. 
This usually implies that we cannot write down the value of Y,,, directly. On the 
other hand a method is called explicit if Y;,, or higher derivatives at x,,, do not 
appear in the recurrence relation for Y,,,. In this case Y,,; is defined explicitly in 
terms of previous values. 


Exercise 1 


Use the trapezoidal method with h = 0.2 to approximate the solution to the differential 
equation 

y’ = 3y + sinx with y(0) = 0, 
at x = 0.2, 0.4 and 0.6. 
The correct solution is 

2s. 1 

yee 10" — 10 

which gives y(0.2) = 0.02460, y(0.4) = 0.12308 and y(0.6) = 0.35304. 


Compare the approximate solution with the true solution by calculating the errors. 


Exercise 2 
Write down the formula to be used when the trapezoidal method with h = 0.1 is to be 
applied to the differential equation 

y’ = siny with y(O) = 1. 
What is the difficulty with this method if you try to use the recurrence relation to generate 
the approximate solution? 
[Solutions on p. 50] 


1.2 Local truncation errors 


So far in this course you have met the Taylor series methods (of which Euler’s 
method is one) and the integration methods (of which Euler’s method is again 
one). By the end of the unit you will be familiar with two more methods. With so 
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many methods available what is the best method to use? This question is by no 
means easy to answer and I shall devote most of Section 3 to examining some of 
the pitfalls. For the moment we will assume that we just have the choice between 
the Taylor series methods and the trapezoidal method and want to make some 
simple comparison between them. 


Ideally for a given problem we would like to compare the error in the nth term in 
the sequence of approximations generated by method A with the error in the same 
term given by method B. 


This is the error 
in the computed 
solution after n steps. 


The point (x,, y,,) has, for 
simplicity, been labelled y, 
and so on. 


Xo) x| X2 X3 poe Xn x 


Figure 3. Graph of the approximate solution and the true solution. 


However, to do this we would need to know the true solution, which is not 
normally available, so we need an alternative method of comparison. There is a 
fairly simple yardstick by which we can judge a numerical method and that is to 
estimate the error introduced in one single step; if we use a method which 
introduces a comparatively small error at each step then the overall error is likely 
to be comparatively small. 


From Unit 2 Section 1 on direction fields we know that there is a whole family of 
solution curves. Thus after r steps we can assume that our approximation Y, lies 

on some solution curve (not necessarily the one defined by the initial condition). 

The error we are going to measure is the distance of Y,, , from this particular 

solution curve (see Figure 4). y 


To illustrate the technique we look at Euler’s method, given by 
Y,41 = 7; + AY; 


Let y = y(x) be the solution curve in Figure 4 on which Y, lies. We have 


(i) é cute OF xX, Xr + x 
(uy) ¥, = m(x,, Y,) = m(x,,Yr) = yr. Figure 4. Graph showing the 
Tin error in a single step. 
¥e+1 =r + hy,. (7) 
The true value, y,,,, can be found by using the Taylor series expansion at x, to 
give 
h2 3 
Veet Ie ee oe tT (8) 
The error in Y,,, is given by Y,,, — y,+, which we can write, using Equations (7) 
and (8), as 
h2 3 
Yaa r+ =O, +) Ty, + hy + ears ioe 
h2 3 
pre fa ee 9 
ae (9) 


This error in Y,,, is called the local truncation error because it is 


(i) local in the sense that we are looking at the error in a single step, i.e. in the 
locality of x,; 
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(ii) a truncation error in the following sense: by comparing Equations (7) and (8) 
it can be seen that we could have obtained the value of Y,,, by truncating the 
Taylor series in Equation (8). 


Look at the local truncation error for Euler’s method (Equation (9)). The first 
term in this series, (—h?/2)y’, is the most important because it contains the lowest 
power of h. For sufficiently small values of h all the other terms become relatively 
small compared with this first term. For this reason we give the first term a special 
name and call it the principal term in the local truncation error. 


When comparing two methods it is the principal terms in their local truncation 
errors which will normally be compared and in particular the power of h in the 
principal term. The method with the higher power of h in its principal term will 
give a smaller local truncation error if h is sufficiently small. 


I must emphasize that the local truncation error does not, in itself, give us any 
idea of the accumulated error after n steps. We shall deduce some results for 
accumulated errors in Section 3. The local truncation error does give us some 
practical information concerning the choice of step-size. Suppose we halve the 
step-size, h, in Euler’s method. Since the principal term involves h’, halving the 
step-size reduces the local truncation error by a factor of 4, approximately. 
Remember, however, that halving the step-size means that we have to do twice as 
many computations to determine the approximation at a given value of x, so that 
there is more accumulation of errors. 


In the following example we shall see how the accumulated error behaves for two 
different values of h. 


Example 2 
Apply Euler’s method to the differential equation 


y=x+y with y(0) = 0, 


using step-sizes h = 0.1 and h = 0.05 to compute the values of Y, at x = 0.1, 0.2 
and 0.3. 


Solution 


0.0025 

0.00763 
0.01551 
0.02628 
0.04010 


The true solution at x = 0.3 is y = 0.0499. The graph in Figure 5 opposite 
indicates the growth of errors. 


Exercise 3 
The Taylor series method of order 2 is given by 
h2 
Text — YT; + hY}. + 7 Yr. 


Calculate the local truncation error for this method. Write down the principal term in the 
local truncation error. 
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0.04 +> 


Pa i 


0.03 


0.02 


0.01 


0 0.1 0.2 0.3 x 


Figure 5. Graph comparing the numerical results for two 
values of h in Example 2. 


Exercise 4 


11 


Calculate the local truncation error for the Taylor series method of order 3. Write down the 


principal term in this error. By what factor would you expect the local truncation error to 


decrease if the step-size were halved? 


Exercise 5 
Calculate the local truncation error for the Taylor series method of order n. 


(The result of this exercise is important and is repeated in the summary.) 


Exercise 6_ 
The Taylor series method of order 2 was applied to the differential equation 


y=xt+y with y(O) = 0. 


Two sets of results were obtained for step sizes h = 0.1 and h = 0.05 for values of x up to 
x = @3: 


0.00125 


0.005127 
0.011764 


0.021025 0.021305 
0.033897 


0.049232 : 0.049696 


The true values are y = 0.005171 at x = 0.1 
y = 0.021403 at x = 0.2 
y = 0.049859 at x = 0.3. 


Compute the errors at x = 0.1, 0.2 and 0.3 for the two step-sizes. By approximately what 
factor have the errors been reduced by halving the step-size? 


[Solutions on p. 51 | 


Although halving the step-size 
reduces the error per step by a 
quarter (approximately), the 
error in the value of Y at 

x = 0.3 has only been 

reduced by half 
(approximately). We 

shall investigate this further in 
Section 3. 
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1.3. Local truncation errors for implicit methods 


Local truncation errors for implicit methods, such as the trapezoidal method, 
cannot be calculated using the techniques you have met in the last subsection. In 
the trapezoidal method we have 


h 
Yar= 2, + ee + Ye+4) (10) 
If we assume that Y, lies on a solution curve we have 
h / / 
Ye+1 =Vr + 5Or + Y+4+1) (11) 


and the difficulty is that we do not know enough about Y;,, to be able to 
compare this formula for Y,,, with the Taylor series expansion for y,+ 1. However 
we can compute the principal term in the local truncation error by replacing Y;., 
by yi, in Equation (11). Although the error in Y;,, affects the other terms in the 
local truncation error, it does not affect the first term since Y;,, is multiplied by h 
in Equation (11). In effect we are saying that the principal term in the local 
truncation error can be computed by assuming that, in Equation (10), all the 
right-hand side values can be replaced by their true values. 


The Taylor series expansion for y+, is 


h2 
Vr+1 = J VIF alee ee 


Substituting this in the approximate expression for Y,,, gives 


h 
t re Pe alr + Yea) 


os = re th oe es (12) 
eas Yr yr 3 r Yr : Yr . 
Now the true value for y,,, is given by the Taylor series at x, as 

2 3 


h 
rst = Ve Ta tee (13) 


Subtracting Equation (13) from Equation (12) gives 


Y. +o eae m \ hy’ h? ” ae 
r+1 Vr+i1—VJ)r avr 7 Yr Vr ar Yr Yr aor 6." ee 
h3 
<< ee 


Thus the principal term in the local truncation error for the trapezoidal method is 
(h°/12)yy". 

The following procedure summarizes the technique for computing the principal 
term in the local truncation error and can be used for both implicit and explicit 
methods. For an explicit method you get all the terms in the local truncation error 
whereas for an implicit method you can only get the principal term. 
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To find the principal term in the local truncation error 


1. Write down the recurrence relation for the method as a 
formula for Y,,,. 


. Replace Y by y in all symbols on the right-hand side of this 
recurrence relation. 


. Expand this new recurrence relation as a power series in h 
using Taylor series about x,. This gives Y,,, in terms of ioe 
eee 

. By subtracting the Taylor series 

2 


/ h “ 
Vr+1 =f 


express Y,,; — y,+1 in powers of h. 


. The term in the lowest power of h is the principal term in the 
local truncation error. 


We can compare, for example, the trapezoidal method with Euler’s method. Since 
the principal terms in the local truncation error involve h? and h? respectively we 
would expect much better results from the trapezoidal method for sufficiently 
small values of h. 


Exercise 7 
Consider a method which uses the recurrence relation 


fee = ¥ + Wl 446 


Determine the principal term in the local truncation error for this implicit method. What 
other method have you seen whose principal term involves the same power of h? 


Exercise 8 

(i) For the trapezoidal method, if you halve the step-size what is the approximate factor 
by which the error per step reduces? Which Taylor series method would reduce the 
error per step by the same factor? 


(ii) Using the results of Exercise 6, by what factor would you expect the accumulated errors 
to reduce if the step-size were halved for the trapezoidal method? 


Exercise 9 
The trapezoidal method was applied to the differential equation 
y=xt+y with y(0) = 0. 


Two step-sizes, h = 0.1 and h = 0.05, were used as in Exercise 6. The following tables of 
results were obtained. 


0.001282 


0.005263 is. 0.005194 
0.011871 

0.021607 ; 0.021454 
0.034092 

0.050197 0.049943 


Comparing these results with the true solution at x = 0.1, 0.2 and 0.3 we have the following 
table of errors. 


; 0.000092 0.000023 
0.2 0.000204 0.000051 
0.3 0.000338 0.000084 


14 MST204 19.2 


By what factor (approximately) have the accumulated errors been reduced by halving the 
step-size? 


[Solutions on pp. 51-52] 


Summary of Section 1 


1. Integration methods for the numerical solution of differential equations are 
derived from numerical integration techniques. 


2. The trapezoidal method is an integration method and gives rise to the 
recurrence relation 


h / / 
Tei = Le a yh = ¥ 544% 


3. The trapezoidal method can be applied to the linear differential equation 


y’ =IU(x)y + k(x) 
using the procedure on page 8. 


4. The trapezoidal method is an example of an implicit method. An implicit 
method for solving a differential equation gives rise to a recurrence relation 
involving Y/,, or higher derivatives at x,,,. A method is explicit if the 
recurrence relation does not involve Y’,, or higher derivatives at x,,,. Implicit 
methods are difficult to use with non-linear equations. 


5. The local truncation error for a numerical method for solving differential 
equations is the error induced in Y,+, assuming that Y, lies on some solution 
curve, y, so that Y, = y,. 


6. The principal term in the local truncation error is the term with the lowest power 
of h. The principal term can be determined using the procedure on page ‘3: 


7. The principal terms in the local truncation error for various methods are as 
follows: 


method principal term 

h2 
(i) Euler’s method sie y" 

h? 
(ii) Taylor series method of order 2 z yn” 

n+1 
(iii) Taylor series method of order n a (n+1) 
: (n+ 1)! 
h 

(iv) Trapezoidal method - y” 


2 The predictor-corrector method 
(Television Section) 


2.1 Pre-broadcast notes 


The television programme for this unit investigates ways of obtaining numerical 
solutions to the logistic equation 


ee i, 
é 101 on 


This equation models the growth in a population where the initial population is 
100 and the population stops growing if it reaches 1000. 


with y(0) = 100. 


This equation has an analytic solution which you will be asked to find in 
Exercise 1, so you may be wondering why we should want to determine numerical 
approximations for this problem. There are several reasons: 


You met this equation in 
Unit 3 Section 2. 
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(i) We can compare our numerical solution with the true solution. 


(ii) If we have to obtain numerical approximations to a closely related 
differential equation which cannot be solved analytically, we may be able to 
use information we derived solving this one. 


(iii) The non-linearity in this equation will highlight some of the difficulties in 
using numerical methods and how these difficulties can be overcome. 


The television programme is divided into three parts: 


(i) The integration methods developed in Subsection 1.1 will be reviewed and 
the methods applied to the logistic equation. An important feature of the 
programme is the way in which the graphs of y’ and y inter-relate. 


(ii) The failure of the trapezoidal method gives rise to a method based on using 
Euler’s method and the trapezoidal method together in a method known as a 
predictor-corrector method. We will look at this method in more detail after 
the programme. 


(iii) The final part of the programme looks at the importance of choosing a small 
enough step-size, h, so that the numerical solution is qualitatively the same 
as the true solution. This topic will be explored more fully in Section 3. 


Exercise 1 
Use the method of separation of variables to solve the differential equation 


y 
Bee eee G 
| i 


What is the particular solution if y(0) = 100? 


Sketch a graph of the particular solution for 0 < x < 1. (The easiest way to do this is to 
plot a few points.) 


[Solution on p. 52] 


Now watch the television programme ‘Integrating by numbers’. 


2.2 Programme notes 


Euler’s method, as an integration method, can be considered as a means of 
building up a graph of the approximate solution, y, by calculating the approximate 
area under the graph of y’, as described in Subsection 1.1. In Subsection 1.1 we 
derived the basic equation for an integration method as 


x 


V1 =Yot | y' dx (from Equation (3) of Subsection 1.1) 


xo 


where y’ = m(x,y) is the differential equation to be solved. 


Euler’s method was derived by approximating the integral by the area of a 
rectangle of width h = x, — Xo and height yo (see Figure 1(a)), 


1.€. | y' dx = hyo. 


The value of yo can be found by substituting the initial condition yo in the 
differential equation to give 
Yo sai m(Xo, Yo): 


The value of the area, hyo is added onto the value of yo to give an approximation 
for y, so that 


Y; = Yo + hyo 
as shown in Figure 1(b). 


Once we have found Y, we can repeat the process to find Y;, then Y, and so on. 
In this way we can build up the graphs of both y and y’. Summarizing the method 
we have a three-stage procedure: 


TV19 


es . aera xX 


Xo a 
Figure 1(a). Graph of y’ 


against x showing the area for 
Euler’s method. 


y 


Yo 
XO xX * 


Figure 1(b). How yo is 
incremented to obtain Y,. 
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(i) Evaluate Y} as Y; = m(x,, Y,) and plot this point on the y’ graph. 
(ii) Calculate the area of the rectangle height Y; and width h = x,,, — X,. 


(iii) Increment Y, by the value of this area to obtain Y,,+;. 


The important thing here is to 
see how the graphs of y and y’ 
inter-relate. 


Stage (111) 


Figure 2. Schematic diagram for Euler’s method. 


If we try to get the same geometrical insight into the trapezoidal method we come 
unstuck because the method requires both y, and Y; in order to be able to 
calculate Y, using the formula 


ae : 
Y; =yot+ 5 (Vo 2 


To find Y, we must increment yy by the area of a trapezium and we cannot find 
the area of the required trapezium geometrically because the only information we 
have is the values of yo and yo = m(Xo, Yo). 


For linear equations we can do some algebraic manipulation to get round this 
difficulty, as described in Subsection 1.1. For the logistic equation we could do 
some algebraic manipulation and determine Y, by solving a quadratic equation. 
This process is much more involved than Euler’s method and in the programme 
this method is not pursued. 


The second part of the programme indicates how we can get round the difficulty 
in using the trapezoidal method for non-linear differential equations by using a 
combination of Euler’s method and the trapezoidal method. To apply the 
trapezoidal method we need Y; in order to compute Y; as 


h 
Y; = Yo + 5 + Y;) 


and Y; depends on Y;. The way round this difficulty is to use Euler’s method to 
obtain a crude approximation for Y; which we can then use in the trapezoidal 
formula. Because this crude approximation is multiplied by h in the formula we 
can still get quite a good approximation, Y;. 


Example 1 


To see how the method works I will compute a value for Y, at x = 0.1 for the 
logistic equation using this new method. 
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Solution 
The method is in four stages: 


(1) 


We know that Yo = 100, so I can evaluate Yj using the 
| Y, ae 
differential equation as Yj = 10Y, f ~ ial giving 
100 
Yo = 10 x 100) 1 — —— 
So | aos 
Yo enables us to calculate the area of the rectangle in Figure 
3(a). 
(ii) | Using Euler’s method I get a crude approximation for Y, in 
Figure 3(b) as 


Y*¥ = % +hY}, = 100+ 01 x 900 
= 190. 


I’ve labelled this with an asterisk to indicate that it is a 
crude approximation. We are going to do better. 
(iii) This is not a particularly good approximation. However I 


can use it to obtain a value for the derivative Y*’, using the 
logistic equation, as 


YY?’ = mn, YT) 
YT 
1000 


190 
Yt = 1o¥t(1 ~ | 


} = 10 x 190{1 — oe 


= 1539. 


This is shown in Figure 3(c). 


Look at Figure 3(d). Since we now have values for Yj and 
Yj’ we can calculate the area of the trapezium. Adding it to 
Yo we obtain an improved approximation Y, as 


(iv) 


h 
eee +l ¥s re rt’ 


0.1 
= 100 + (900 + 1539) 


= 721.95, 
as shown in Figure 3(e). 
This is much closer than Y# to the true value of y at 


x = 0.1, which is 231.97 correct to 2 decimal places. 


The method outlined above is an example of a predictor-corrector method. In 
Example 1 Euler’s method was used to predict a value for Y,. This value for y, 
was then corrected using the trapezoidal method. 


Predictor-corrector methods provide one of the most widely used techniques for 
solving differential equations. They are based on the idea of using an explicit 
method to obtain a prediction and then using a more accurate implicit method to 
correct this prediction. 


The four stages of the procedure I used in Example 1 can be summarized as: 
(i) 

(il) 
(i11) 
(iv) 


Evaluate Yo using the differential equation. 
Predict Y¥ using Euler’s method. 
Evaluate Y*' using the differential equation. 


Correct Y, using the trapezoidal method. 


l7 


x; = 0 x,=01 x 


Figure 3(a) 


Figure 3(c) 


Y : , 
Yi 


5 (Yo + Y7) 


0 0.1 
Figure 3(d) 


0 (0.1 x 
Figure 3(e) 
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Here is the same procedure in diagramatic form. 


(iii) EVALUATE 


(i) EVALUATE 


Area is 


5 (Yo'+ 71) 


(iv) CORRECT 
y 


TRAPEZOIDAL 
Y= Yo+5 (Y5+¥1) 


EULER 
Y; = Yo a hY} 


Figure 4. Schematic diagram of the Euler-trapezoidal predictor-corrector method. 


Example 2 

When the procedure above was used to compute Y, for the logistic equation, the 
value Y, = 221.95 was obtained. I can now use this value as the starting point in 
the computation of Y>. 


Solution 
We can use the four stages again as follows: 


(i) Evaluate Y; = 10, (1 — Sie) = 10 x 221.95(1 — — 
= 1726882, 
(ii) Predict Y$ using Euler’s method as 
YS =:¥, + hY; 
= 221.95 + 0.1 x 1726.882 
= 394.6382. 
(iii) Evaluate Y¥’ = 10Y3 [i — iio) = JO x 394.6382 1 — ae | 


= 2388.9889. 


(iv) Correct Y, using the trapezoidal method. 
ae 
Y, = Y; BG + Y#’) 


= 221.95 + 0.05(1726.882 + 2388.9889) 
= 427.74355. 


Exercise 2 
Follow the method of Example 2 to compute values for Y; and Y, for the logistic equation. 


[Solution on p. 52] 


Thus we can use the Euler-trapezoidal predictor-corrector method to compute ¥,, 
Y>,..., Yq. aS summarized in the following procedure. 


hY; is the area of the 
rectangle. 


h 
a (Ys + Y#’) is the area of the 


trapezium. 
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The predictor-corrector method (Euler-trapezoidal) 


1. You are given a differential equation 
y’ = m(x, y) 
and a value for y(Xo). 


. To apply the predictor-corrector method, choose a step-size h, 
and calculate Y,, Y,, Y3,... by repeating the following cycle of 
steps starting with Yo = y(Xxo): 


(i) Evaluate Y, = m(x,, Y,). 

(ii) Predict using Euler’s method: 
Y*.1,= Y, + hyY,. 

(iii) Evaluate Y*4,, = m(x,, Y*,;). 


(iv) Correct using the trapezoidal method: 


h 
Yen = Vetegtle® Ue). 


3. Y, is an approximation to y,. 


Note that any predictor-corrector method is an explicit method. (You will be 
asked to show this for the Euler-trapezoidal method in Exercise 5.) 


In general, predictor-corrector methods use an explicit method to predict Y*?, , 
and a more accurate implicit method to correct Y,4,. 


Although the derivation of the local truncation error for a predictor-corrector 
method is beyond the scope of this course, it can be shown that the principal term 
in the local truncation error for the Euler-trapezoidal method involves h* (as does 
the principal term in the local truncation error for the trapezoidal method). 


The following table gives the computed values for the logistic equation for each of 
the three methods discussed in the programme and also the true solution, obtained 
using the formula from Exercise 1. The results for the trapezoidal method were 
obtained by solving a quadratic equation at each step. 


Solutions to the logistic equation, with h = 0.1. 


trapezoidal predictor-corrector true solution 
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I have plotted these points on a graph shown in Figure 5. Note that the 


trapezoidal method gives the best results with the predictor-corrector method 


initially giving much better results than Euler’s method. 


ion 


ic equat 


1st 


Figure 5. Comparison of results for the log 
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In the final part of the programme we look at what happens if we use different 
step-sizes, h. The first method we look at is Euler’s method. The four graphs in 
Figure 6 show the results for h = 0.1, 0.15, 0.2 and 0.25. The logistic curve is also 
shown for comparison. 


y y 
1000 1000 
800 800 
600 600 
400 400 

200 

0 02 A. “26 -O8 = 10 ¢ 0 

y y 
1000 1000 
800 800 
600 600 
400 400 
200 200 i 

0 0.4 0.8 1.2 1.6 2.0 x 0 


Figure 6. Euler’s method for different step-sizes. 


The results show that, for h = 0.1 and 0.15, the approximate solution at least has 
approximately the same qualitative behaviour as the analytical solution. The saw- 
tooth behaviour near y = 1000 of the solution for h = 0.2 and 0.25 is clearly 
undesirable. This will be discussed in Section 3 but the graphs do indicate that, for 
some problems at least, we have to be careful about our choice of step-size if we 
want to get reasonable results. 
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The trapezoidal method, on the other hand, behaves well for this problem as the 
graphs in Figure 7, with step-sizes h = 0.1, 0.2, 0.5 and 1.0, illustrate. 


1000 
800 
600 
400 


200 


1000 
800 
600 
400 


200 


0 0.5 1.0 3 2.0 6 He SE 2 0 


Figure 7. The trapezoidal method for different step-sizes. 


While numerical results for h = 0.5 and 1.0 are not very good they do converge to 
the limit y = 1000 as x increases and reasonably approximate the shape of the 
logistic curve. 


So what about the predictor-corrector method? It is a combination of the Euler 
and trapezoidal methods. Unfortunately, as indicated in the graphs in Figure 8, 
this method also behaves badly for larger values of h. 
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600 600 


400 400 
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1000 1000 

800 800 

600 600 

400 400 
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Figure 8. The predictor-corrector method for different step-sizes. 


The numerical results follow the shape of the logistic curve for h = 0.1 and 

h = 0.15, but for h = 0.2 the numerical approximations are converging extremely 
slowly towards the limit value of 1000 while for h = 0.25 the numerical results 
converge rapidly to the wrong limit. 


Exercise 3 


Use the predictor-corrector method (Euler-trapezoidal) with h = 0.1 to determine 
approximate solutions for y at x = 0.1 and x = 0.2 for the differential equation 


y =logy+x with y(O) = 1. 
Quote your results to 3 decimal places. 


Exercise 4 


Use the predictor-corrector method (Euler-trapezoidal) with h = 0.2 to determine 
approximate solutions for y at x = 0.2 and x = 0.4 for the differential equation 


y’ = siny with y(O) = 1. 
Quote your results to 3 decimal places. 
Exercise 5 
Eliminate Y*,, and Y*;, from the recurrence relations for the predictor-corrector method 


P54 = FF BY, where Y= m(x,, Y,) 
h 
and Kai = Yet git + Fea where YF41 = m(x,41, Y7+1) 


and hence verify that this method is an explicit method. 
[Solutions on pp. 52-53 ] 
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Summary of Section 2 


1. The television programme looks at several methods of obtaining numerical 
solutions to the logistic equation 


1000 
including Euler’s method and the trapezoidal method. 


y’ = 10y (1 a with y(0) = 100 


2. A predictor-corrector method uses an explicit method to obtain a crude 
approximation to the solution of a differential equation at x,,, and then uses a 
more accurate implicit method to obtain an improved approximation at x,+1. 


3. The predictor-corrector method based on Euler’s method and the trapezoidal 
method overcomes the difficulty in using the trapezoidal method to obtain 
numerical solutions for non-linear differential equations. The procedure for 
applying this method is given on page 19. The principal term in the local 
truncation error for this method involves h°. 


4. For Euler’s method and the predictor-corrector method we need to choose the 
step-size, h, carefully in order to obtain reasonable results. 


3. The analysis of numerical methods 


From the empirical results of Sections 1 and 2 we have seen that 


(i) it is possible to approximate differential equations using a variety of different 
methods, each of which leads to a recurrence relation; 


(ii) the smaller the step-size, h, the better the approximation; 
(iii) larger values of h may give nonsense. 


In this section we want to look at properties of methods theoretically so that you 
can evaluate any new method you might meet. The analysis will also help you to 
choose the step-size required to give reasonable results. 


We shall introduce three new concepts: 


(1) consistency, which is related to (i) above, 
(2) convergence, which is related to (ii) above, 
(3) stability, which is related to (iii) above. 


For simplicity we are only going to consider one-step methods in which the 
differential equation is approximated by a first-order recurrence relation. All the 
methods you have met so far have been one-step methods, but in Section 4 you 
will see a two-step method in which the recurrence relation is of second order (i.e. 
the formula for Y,,, involves Y,_, as well as Y,). The analysis can be extended to 
multi-step methods but we do not have time to do so in this unit. 


3.1 Consistency 


One of the fundamental requirements of any numerical method that we use to find 
approximate solutions to a differential equation is that the recurrence relation to 
be used should be a reasonable approximation to the differential equation itself. In 
other words the recurrence relation must be consistent with the differential 
equation. If this requirement is not satisfied we cannot hope that the solution 
obtained using the recurrence relation is a reasonable approximation to the 
solution of the differential equation. 


As in the definition of the local truncation error, consistency is concerned with the 
local behaviour of the method. We assume that Y, lies on some solution curve so 
that Y, = y,. For implicit methods we also assume that any other values on the 
right-hand side of the recurrence relation are correct. 
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Any one-step method, including the predictor-corrector method of Section 2, must 
involve a recurrence relation of the form 


oe = is + he(x,, i. Y41;2) (1) 


where the function @ depends on the method used. For example the Taylor series 
method of order 2 uses the recurrence relation 


h2 
i 3 — ee Te 


h 
=. +afr, et 
and so, for this method, we have 
oe 
(x,, te eee h) = i, + 2 pat 


Using the condition that Y, = y, and that the other right-hand side values are 
correct in Equation (1), we have 


ee = Ss = NO Ve Vr+4: h). 
Rearranging this equation gives 
D4 1 


= 1X Vrs Vy 41,10. 3 (2) 


The left-hand side of Equation (2) represents the slope of the line joining the 
points (x,,y,) and (x,41, Y,41). As h gets smaller and smaller we require that 
Equation (2) should be a better and better approximation to the differential 
equation 

y' = m(x, y) (3) 
at x,. This condition holds when h = 0 provided that 


P(X Vr> Ver, 0) = m(x,, y,)- 


(Since y,,, = y(x, +h) we have y,.,; = y, when h = 0.) We thus have a fairly 
simple test for consistency. 


A one-step method is said to be consistent with the differential The definition of consistency 
equation (3) if the recurrence relation can be expressed in the can alternatively be stated as 
ti 1 i h = i, see S 
form of Equation (1) wit lim = Vr oo 
h>0 


P(X Vis Hp, OF = ey): 


Example 1 
For the Taylor series method of order 2 we have 


h 
P(X Vrs Vea 12) = Vy, Be ar 


so that putting h = 0 in this equation gives 
P(X+,Vr, Vr, 9) =), 
= m(Xx,, Yr) 
as required. 


Hence the method is consistent with the differential equation (3). 
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Exercise 1 
Which of the following methods are consistent with the differential equation 


y' =m(x,y) ? 


h 

(i) V+ ee Y, a qtr + a2 5% 1) 
ee h f 

(ii) Y,+1 ae q, 3 Pea is Y ++ 1). 


Exercise 2 
What condition must be imposed on the values of a and b for the method given by 


Y,41 = Y, + h(aY;, + bY,+1) 
to be consistent? 


[Solutions on p. 53 | 


3.2 Convergence of one-step methods 


In Exercises 6 and 9 of Section 1 you were asked to compare the results of using a 
method with two different step-sizes. In Exercise 6 you saw that halving the step- 
size for the Taylor series method of order 2 reduced the errors by a factor of 
approximately 4. In Exercise 9 you obtained a similar result for the trapezoidal 
method. 


What would have happened if we had halved the step-size again? Well, 
commonsense tells us we could expect a further reduction in the error by another 
factor of 4 each time. Thus if we continue halving the step-size (and doubling the 
number of calculations) we would theoretically obtain answers to any accuracy we 
require. (In practice there is a limit to the accuracy that we can achieve because of 
the build-up of rounding errors in the calculations.) Suppose we assume that all 
calculations are carried out exactly. Then our empirical evidence from Exercises 6 
and 9 in Section 1 indicates that the solutions converge in each case. That is, if we 
take some particular point x* (say x* = 0.2) and compute the approximation to y 
at that point using smaller and smaller step-sizes, then these approximations get 
closer and closer to the true solution at that point. Figure 1 illustrates how the 
numerical solution might converge to the true solution. 


y(0.2) 


-—_- 
-—_ 


ad 
- 
-—_ 


Figure 1. Graph showing how the numerical solution could converge to 
the solution of the differential equation. 


For a given numerical method let us look at the error at a fixed value of x as h 
changes. 


(i) Suppose we take a fixed value of x, say x*. 
(ii) As h decreases so the number of steps, N*, from Xo to x* increases such that 
x* = Xq + N*h. 
The global error at x* is the error 
Yv«e — y(x*). 


We are now in a position to define convergence. 
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1. You are given a problem of the form 


y = m(x,y) with y(x9) given. 


2. You wish to use a one-step numerical method to obtain 
approximations for x») <x <b. 


3. The method is convergent for this problem if, for all points x* such 
that xp < x* < b, we have 


lim Yy+ = y(x*), where N*h = x* — Xo 
h-0 


i.e. the global error at x* tends to zero as h tends to zero. 


Convergence is difficult to prove directly, except for very simple differential 
equations. However we can state a rather simple but very useful theorem for one- 
step methods which makes any direct proof of convergence unnecessary. 


Theorem 1 
For a given problem 
y’ = m(x, y) 


and a value for y(x 9), where m is well-behaved, a one-step method 
is convergent if and only if it is consistent. Furthermore if the 


principal term in the local truncation error involves h?**, for 
some integer p, then the global error at x*, for small h, is of the 
form 


Yy« — y(x*) = Ch? 


where C does not depend on h. 


This is an important theorem because it gives us a way ot evaluating our methods. 
For example, the trapezoidal method has a local truncation error whose principal 
term involves h°. Our theorem tells us that the global error is of the form 


Yn« = y(x*) = Ch. 


While this in itself does not help to put a bound on the error (since we do not 
know C) it indicates that if we were to halve the step-size h then we could expect 
the error in Yy+ to decrease by a factor of four, even though we’ve had to do twice 
as many calculations to compute Yy«. This is precisely the result we deduced in 
Exercises 6 and 9 in Section 1. 


The above results indicate that we could classify methods in terms of the principal 
term in their local truncation errors. A preliminary rank ordering of the methods 
we have met so far, indicating their practical limitations, would be as follows. 


Well-behaved here means that 
we can differentiate m as often 
as we wish. 
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Principal term 
involves 


Taylor series method 
of order 3 
Drawbacks: 

second and third 
derivatives 

required. 


Predictor-corrector Taylor series method Trapezoidal method 
(Euler-trapezoidal) of order 2 


Drawbacks: Drawbacks: Drawbacks: 
none second derivative for linear 
required. equations only 


Euler’s method 


Drawbacks: 
none. 


This diagram makes the point that if high accuracy is required then the best 
method is likely to be a Taylor series method provided the derivatives can be 
obtained. If less accuracy is required or the derivatives are difficult to compute 
then the predictor-corrector method may be preferred. 


Exercise 3 


Show that the trapezoidal method is consistent and hence convergent. How does the global 
error behave at a point x* for small values of h? If the step-size were halved, by what factor 
would the global error at x* be reduced? 


Exercise 4 


Use the results of Exercise 5 in Section 2 to show that the predictor-corrector method 
(Euler-trapezoidal) is consistent. 


[Solutions on p. 53 | 


3.3 Stability 


There is still an outstanding practical issue in the analysis of one-step methods in 
that we do not yet know, given a particular differential equation, what value for 
the step-size, h, to use. The convergence theorem has only told us that, as h tends 
to zero, the numerical solution tends to the true solution—and only then if we 
carry out our computations accurately. For a chosen step-size we need to know 
whether we can expect reasonable results. 


In the television programme we saw examples of how numerical solutions can go 
haywire if the chosen step-size is too large. In each case problems arose because 
errors, incurred in approximating the differential equation, swamped the solution. 
In this subsection we are going to investigate a fairly simple test equation e see 
when these difficulties arise and in the next subsection we will apply our results to 
more practical problems. The following exercise illustrates what can happen. 


Exercise 5 
Use Euler’s method with h = 0.1 to solve the differential equation 
y’ = —30y _ with y(O) = 1. (1) 


Compare your results with the true solution, y = e~ °°. 


[Solution on p. 53 | 
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In Unit 1 two types of ill-conditioning were defined: absolute ill-conditioning and 
relative ill-conditioning. In this unit we are only going to discuss absolute ill- 
conditioning. 


First we consider the differential equation itself. Some differential equations are 
absolutely ill-conditioned in the sense that small changes in the data result in 
much larger changes in the solution. Whatever numerical method we use we 
cannot hope to avoid the build-up of errors, although we can hope that those 
errors will not swamp the solution. In order to deal‘with this case thoroughly we 
would need to consider relative ill-conditioning, which is beyond the scope of this 
unit. On the other hand there are some differential equations which are well- 
conditioned and it is these problems which particularly concern us here. 


When we apply a numerical method to a well-conditioned problem the recurrence 
relation problem can be absolutely well- or ill-conditioned. In order to discuss 
this we need a new definition which applies to numerical methods rather than to 
the problems they are intended to solve. 


A numerical method, applied to a differential equation, is 
absolutely unstable if the resulting recurrence relation problem is 
absolutely ill-conditioned. The method is absolutely stable if the 
recurrence relation problem is absolutely well-conditioned. 


To illustrate these ideas we consider the test problem given by 
y =ay with y(0) = 1 


where « is a known constant. This equation can be solved analytically (using, for 
example, the separation of variables method) giving the solution 


y=e" 


Figure 2 shows the solution for various values of «. 


y 
a2 
o=] 
a=0 
a= —] 
a= -—2 
0 0.5 1 x 


Figure 2. Graphs of e** for various values of «. 


Suppose we perturb the initial condition for this differential equation so that 
y(O)=1+e 

where ¢ is small. The new solution is 
y =(1 + eje* = e* + ce™ = y + ee™ 


The error in the solution (ge**) is significantly larger than the error in the data (é) 
whenever e** is significantly greater than one. In this case the differential equation 
‘is absolutely ill-conditioned. Now for sufficiently large positive x this condition 
holds whenever ais positive. We conclude that the differential equation problem is 
absolutely ill-conditioned if « > 0. By a similar argument the problem is absolutely 
well-conditioned if « < 0. 
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So what happens if we try to use a numerical method to solve the test problem, 
Equation (1)? To illustrate the approach we look at Euler’s method which uses the 
recurrence relation 


V04 = 7, +91, 
From the differential equation, Y; = «Y,, so we have 
Vy = 7, + hey, 
= (1 +ha)Y,. 
We have a recurrence relation of the form 
Y,+1 = @Y, (2) 
where a = 1 + ha. 


The recurrence relation problem is thus to generate a sequence of numbers 
Y,, Y>, Y3... using the recurrence relation (2) with Yo = 1. 


What we now ask is ‘when is this problem absolutely ill-conditioned’?” 
From Unit 1 we know that absolute errors grow whenever 

igi > 1. 
(For example if we put Yo = 1 + « then 

Y, = (1 + s)a" = a" + ea" = Y, + €a" 


and the absolute error, ¢, in Yo has induced an error in Y, of ea". This error is 
larger than the error in the data if |a| > 1.) 


There are two cases to consider: 
Case (i) « > 0 


Since the step-size h is positive we automatically have 


a=1+ha>1 y 


true 
and the recurrence relation problem is absolutely ill-conditioned. This is not solution 
surprising as the original differential equation problem is absolutely ill- 

conditioned. However, this ill-conditioning may not matter as the growth of the 


relative errors may be tolerable. 


Case (ii) « < 0 


This is the more interesting case since we know that the differential equation is Figure 3. While absolute _ 
absolutely well-conditioned errors are increasing, relative 
: errors may not increase 
For negative « we have significantly. 
a=1+ha <1. 


However this does not imply that the method is absolutely stable. For example, 
suppose that h and « are such that ha = — 3. Then 


a=1+ha=1-3= -2 
and, since |a| > 1, the recurrence relation problem is absolutely ill-conditioned. 
In general the condition for absolute ill-conditioning is that the modulus 

11 + ha| > 1 


so that whenever 1 +ha < —1 we have an absolutely ill-conditioned recurrence 
relation problem. The condition for absolute ill-conditioning is that 


ha < —2. 


Here we have an example in which the original problem is well-conditioned but 
the numerical method gives rise to a problem which is absolutely ill-conditioned if 
ha < —2. 
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We conclude from this test problem that Euler’s method is absolutely stable if h 
and « are such that 


~2<ha <0. (3) 


The set of all values of ha which satisfy this inequality is denoted by (—2,0) and is 
called the interval of absolute stability for Euler’s method. 


We can now use this analysis to understand what went wrong in Exercise 5. There 
we tried to solve the differential equation 


y' = —30y with y(0) = 1 


using Euler’s method with a step-size h = 0.1. The theory tells us that the method 
is absolutely stable for this problem only if 


—2 < —30h < 0, 


1.e. h < 1/15. While the theory does not tell us that we will get accurate answers, 
for h = 0.05 say, it does warn us that it is impossible to even get sensible answers 
for the step-size we tried, h = 0.1. 


The tape session in the next section will help you to prove some of the stability 
results for other methods. 


Exercise 6 
For what step-sizes is Euler’s method absolutely stable for the differential equation 
y’ = —100y with y(0) = 1? 
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3.4 Stability in practice 


I mentioned in the last subsection that the absolute stability results we obtained 
for the equation 


y =ay with y(0) = 1 


could be extended to more difficult problems. Indeed the extension to the general 
linear differential equation 


y’ = I(x)y + k(x) 
is quite straightforward. Euler’s method gives the recurrence relation 
Y,41 = Y, + hY, 
= Y, + h[l(x,)Y, + k(x,)] 
= [1 + Al(x,)]Y, + hk(x,). 


This non-constant-coefficient recurrence relation is absolutely stable (see Unit 1, 
Subsection 3.2) if 


|1 + Al(x,)| < 1, 
1.e. —2 < hl(x,) < 0 for all the values of r. 


We would thus have to check, before implementing the method, that this 
relationship holds for all the values of x, we are going to use. However, if I(x,) > 0 
we can still use Euler’s method with some confidence provided that we know that 
the solution is increasing in magnitude. 


The treatment of non-linear equations is not quite so straightforward. We first 
consider equations of the form 


y' =m(y). 


For example 


y = 
oT with y(0) = 100. 


(This is the problem considered in the television programme.) 


y = 1051 


In M101, Block 1, the 
notation ja,b[ was used for 
the open interval a < x < b. 
Here (— 2,0) means the same 
as ]—2,0[. 


This is the same as Equation 
(3) with « replaced by /(x,). 
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If we apply Euler’s method to this type of problem we obtain the recurrence 
relation 


Y3 = T, a hm(Y,). (4) 


To see whether errors grow, we look at what happens to Y,,, if there is an error, 
é, in Y,, so that the computed value for Y, is 


Y= ¥,+<6. 
Then the computed value for Y,,, is 
Y¥41= ¥+hm(Y) 
= Y,+¢é+ hm(Y, + 8). (5) 


We can approximate m(y) near Y, by its Taylor polynomial of degree one as 


mY, + 6) =~ m(Y,) + heey 
dy 


Hence, using this approximation in Equation (5), we have 


dm 


Yu, Y,+e+hm(Y,) + he 
dy 


(Y,). 


Using Equation (4) we can reduce this to 


dm 
a ie 1 = her}. 
dy 
So the error in the computation of Y.,, is approximately 
d 
1 + non (Y,)}. 
dy 
This means that the absolute error, ¢, in Y, has been magnified by a factor 


f + _ 3) in the computation of ¥.,,. The recurrence relation problem for 
y 


Euler’s method is thus absolutely ill-conditioned if 


4, 


is fats 
dy 


Hence Euler’s method is absolutely stable if 
ee (6) 
dy 


If this condition is true for all values of Y, in the calculation then Euler’s method 
will be absolutely stable for the differential equation y’ = m(y). 


Example 2 
For the logistic equation considered in the television programme we have 
mis) @ 1 = 
1000 
Differentiating with respect to y gives 
d 
= ii 
dy 50 


d 
Now y takes values between 100 and 1000 so that a takes values between 8 and 
y 
— 10. 
Our conclusions would be as follows: 


; dm 
(i) For 100 < Y, < 500, a 
y 
_unstable for any value of h. However, since the solution is increasing rapidly in 
magnitude, this instability may not cause significant errors in the solution. 


(Y,) is positive and so Euler’s method is absolutely 


This is the same as Equation 


d 
(3) with « replaced by =(%) 
y 
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d 
— (Y,) is negative and can be as small as — 10. 


(1) For 500 < Y, < 1000, 
dy 


d 
Equation (6) is satisfied with = = —10 only if h < 0.2. Thus we would only 
y 
contemplate using Euler’s method with step-size h < 0.2. 


The above analysis agrees precisely with the results we obtained in the television 
programme. The instability is most serious for values of Y, near 1000. In Figure 6 
of Section 2 we saw that, for values of h less than 0.2, the sequences generated by 
the recurrence relation converged to 1000 while, for h > 0.2, the sequences 
generated did not converge to 1000. This was caused by the absolute instability of 
Euler’s method for h = 0.2. 


What we have found so far, for non-linear equations of the form 


y =m(y), (7) 
is that the stability condition on the values for h in Euler’s method can be 
——— d 
obtained by substituting the most negative value of - 
y 


—2<ha <0. 


for a in the inequality 


It can be shown that, for any numerical method applied to the differential 
equation (7), the stability condition on the values of h can be obtained by 


substituting the most negative value of — for « in the inequality which determines 
the interval of absolute stability for that method. 


Finally we consider the general equation For example 


y’ = m(x, y) t= 


with y> given with (0) = 1 for 0 < x < 10. 
0 : 


If we apply Euler’s method to this problem we obtain the recurrence relation 
Y-+1 = Y, + hm(x,, Y,). 
This can be written as 


Drei = Y; + hm,(Y,) (8) 


: ; F 2 eee : tx,= 3 
where m, is the function defined by oS Moe red 3y2 = a 


function of y only. 
m,(y) = m(x,,y). 


The treatment of Equation (8) is now identical to that given to Equation (4) and 
we can deduce that Euler’s method is stable if 


dm, 
dy 


for all values of r. 


—2<h 


(77 =0 


Hence, when we consider the stability of a numerical method applied to the 
differential equation 


y’ ie m(x, y), 
we replace « in the inequality which determines the interval of absolute stability by 
dm 
dy 
previous examples are special cases of this general first-order differential equation.) 


“(Y,) where m,(y) = m(x,,y). (Note that all the 


the most negative value of 


Example 3 
For the equation 


y’ = —xy? with y(0) = 1 
for 0 < x < 10 we have 


m,(y) = —x,y” 
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dm 


“~(Y,) = —2x,Y,. 
7 id = 28 


and 


From, for example, a sketch of the direction field for this equation we could show 
that y lies between 0 and 1. Thus, since 0 < x < 10, the most negative value of 


—2x,Y, is —20 and Euler’s method is certainly absolutely stable if If we knew the solution we 
could get a less restrictive 
—2 < —20h <0 bound on h but we have 
; managed to obtain a useful 
ie. h <2, bound on h even without 


knowing the solution. 
From a practical viewpoint, faced with a differential equation, a simple procedure 


for obtaining numerical solutions to a required accuracy is to 


(i) determine bounds, if any, on the values of h to ensure the absolute stability of 
the method; 


(ii) solve the differential equation problem using the numerical method with two 
different step-sizes, h. A comparison of the results at specific x values will 
normally indicate the accuracy. If the discrepancies are larger than desired, a 
smaller step-size should be used and the numerical results compared. 


Exercise 7 
Euler’s method is to be used to solve the differential equation 


1 
y ai eae ee with y(0) = 1 


for0 <x < 10. 


Determine a bound on the values of h which will guarantee that Euler’s method is 
absolutely stable. 


Exercise 8 
Euler’s method is to be used to solve the logistic equation 


— 2. 
1000 
with the initial condition y(0) = 2000. 


y = 105 


The theory tells us that the solution should tend to the stable population of 1000. 
Determine a bound on the values of h which will guarantee that Euler’s method is 
absolutely stable. 


Exercise 9 
Euler’s method is to be used to solve the differential equation 


6 
fees with y(0) = 1. 
y 


It is known that y takes values between 1 and 5 for 0 < x < 2. Determine a step-size h for 
which Euler’s method will be absolutely stable. 
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Summary of Section 3 

1. A one-step numerical method, which uses the recurrence relation 
Ya3 = % +066, Ta hea) 

is consistent with the differential equation 

y’ = m(x, y) 

if 
P(Xrs Vrs Vrs O) = M(X,; Yr). 

2. The global error at x* is the error 
Yy« — y(x*) 


where y(x*) is the true solution at x* and Yys is computed using a recurrence 
relation N* times, where 


: Xo + N*h. 
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3. A one-step method is said to be convergent for the problem 
y’ = m(x, y) with y(x9) given, 


if the global error at x* tends to zero as the step-size, h, tends to zero, for all 
points x* in the interval, [X9,5], in which approximations are required; 


1.€. lim Yy + = y(x*), 
h>0 


where N*h = x* — Xo. 


4. A one-step method applied to a well-behaved differential equation is convergent 
if it is consistent (see Theorem 1). If the local truncation error involves h?*! 
then the global error at x*, for small h, is of the form 


Yy « — y(x*) = Ch? 
where C does not depend on h. 


5. A one-step method, applied to a differential equation, is absolutely unstable if 
the resulting recurrence relation is absolutely ill-conditioned. The method is 
absolutely stable if the recurrence relation problem is absolutely well- 
conditioned. 


6. The interval of absolute stability (a,b) is the set of all values of ha for which the 
method is absolutely stable. 


7. For Euler’s method we have the following intervals of absolute stability 


differential equation | interval of absolute stability 


y’ =ay —2<h«<0 
y’ = I(x)y + k(x) —2<hli(x,)<0O for allr 
y =my) —2< noe) =< foralir 
dm, 
y’ = m(x, y) —2<h—(¥,)<0 — where m,(y) = m(x,, y) 


dy 
For other methods it is sufficient to determine the interval of absolute stability for 
the test problem 

y= oy with y(0) = 1 


and then to deduce the appropriate conditions for other differential equations as 
above. 


4 Exercises on stability and Simpson’s 
method 


4.1 Exercises on stability (Tape Subsection) 


In Subsections 3.3 and 3.4 we looked at the absolute stability conditions for 
Euler’s method. In the first part of the tape session we see how the intervals of 
absolute stability are determined for the trapezoidal method and the predictor- 
corrector method. 


In the second part of the tape we apply these results to the logistic equation of 
Section 2 to confirm the empirical evidence found in the television programme. 


Start the tape when you are ready 


x 
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Recurrence relation for the trapezoidal method 


herve 
Mere x te a cy, to) 


lal <1 when he < aS 


The interval of absolute stability for the trapezoidal method is (—00, ow 
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yo =“y y(0O) =1 


Predictor - corrector method (Euler - trapezoidal) 


trapezoidal 


4 See eee ae 
Pee ge Se | oe a 


case saccsescaeuee snes “ted sscccessen vanne sane Bue tet aes = oa 


} pasense Seanene 
speutaabedtatad ates afentanteteatessettitattateats 
SEES EEE HaGE eatHaEE Hn SHES tta 


saeiesisieaitissiis cit HEHE HHH 


Spassaceataaescaesersececerser caceaeneereae 

Sei ee eg 

ae itt 

Geeeeee eer oe SESE EERE fe HEE 
saat ar ear ay ay eterna Hee 
feritteaieeceestttst ea oa ues Eaiaia aie eta 
eescesegnt doeeneeseces motte th SeSSReRaRGREEESSt seseraseaeeeesss CCC 


is ai aaes S55 8a see seas eecaai SESE iiasseeeeescees, 

HEE HERE, THREE, pease 5 ies aera ta 
Bienats Seunueceess SSUCRRSGERREEEGEAS Feeeue seieadttoitiast ee 

siseriseesee Hs FH aHRECadi oaaiasai fossa ieeaiitaatittaatt seerratirini: HH sagusetuscessesieris 


ranean! 


RRSRS SASSER SS e 


la|<| when| |< he<| | (4a) 


The interval of absolute stability for the predictor - corrector 


method (Euler - trapezoidal ) is a 
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The logistic equation 


y= 10y ee y (0) = 100 


If 500 < Tis 1000 the minimum value of “y is =e 


The trapezoidal method 


The interval of absolute stability is (—0, 0) 
. the method is absolutely stable if he < 0. 


For the logistic equation, with 500 <y < 1000, what 


values of h make the trapezoidal method abeolutely stable ? 


cae 


The predictor~corrector method (Euler- trapezoidal) 


The interval of absolute stability is (-2, 0) 
The method is absolutely stable if -2 < ha< 0 


For the logistic equation the method is absolutely stable 
for 500 <y < 1000 if 


dm 
toms ——  < 
eae 0 


~2 = 10h 
ie h <0-2 
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Exercise 1 
Find the interval of absolute stability for the Taylor series method of order 2, applied to the 
test equation y’ = ay with y(0) = 1, given by 

2 


h 
Yi1=¥,+hY+>¥7 


where Y” = aY) =a7Y,. 
For what values of h is this method absolutely stable when applied to the logistic equation 


with y(0) = 500? 


y 
ewe 
3 | ma 
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4.2 Simpson’s method 


For simplicity we have only analysed one-step methods so far. Most of the 
methods used in practice are many-step methods. To give you some idea of how 
they work we are going to consider a two-step method: Simpson’s method. Like 
the trapezoidal method Simpson’s method is an integration method. Because of its 
high accuracy Simpson’s method was very widely used, before the advent of 
computers, to solve differential equations. It is possible to show that the local 
truncation error for Simpson’s method involves h° so that the global error using 
this method is likely to be much smaller than for most of the methods we’ve 
already discussed. 


Simpson’s method approximates an integral over two steps as 
Xr+i1 h 
f(x) dx = 5 (f(%r-1) + 4f (x,) + f(x,+1)). (1) 
We are going to use this formula to solve the differential equation 
y’ = m(x, y). (2) 
If we integrate both sides of Equation (2) over the interval [x,_,,x,+1] we have 


[yay= | monyyax 


Xr-1 Xr-1 
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The left-hand side is just y,4; — y,—1 So that, using Simpson’s method (1) for the 
integral on the right-hand side, we have 


yeni Pent | m(x, y) dx 


h 
= E (m(x,- 1>Yr- 1) 2 4m(x,, Yr) — m(X,+ 1>Yr+ 1)). 
This leads us to define the recurrence relation for Simpson’s method as 
h 
Y,41=%-1 + qtr i AT 2) 


where Y;, = m(x,, Y,). 


This is a second-order recurrence relation (because it involves terms at x, and 
x,—,1) and it is also an implicit method (because it involves terms in Y’,, on the 
right-hand side). This causes two difficulties: 


(i) Second-order recurrence relations require two initial values, Yo and Y;, in 
order to generate the sequence Y>, Y3, Y,... . However, we only have one 
given initial condition, Y), so that somehow we have to generate a value for Y, 
in order to use Simpson’s method. This is done by using another method, such 
as a Taylor series method, to generate Y, before switching to Simpson’s 
method. 


(ii) Simpson’s method cannot be easily used on its own to solve non-linear 
equations because it is implicit (cf. the trapezoidal method). For this reason 
Simpson’s method is more commonly used as a corrector in a predictor- 
corrector method. 


By a process of algebraic manipulation it is possible to obtain a linear recurrence 
relation for Simpson’s method applied to the linear differential equation 


y’ =I(x)y + k(x) 
with y(xo) given. 


I will not present those manipulations here but the process is similar to the 
procedure for the trapezoidal method in Section 1. The result is given in the 
following procedure box. 


Simpson’s method for linear differential equations 
1. You are given a linear differential equation 
y’ = I(x)y + k(x) 


and a value for y(Xo). 


2. To apply Simpson’s method, choose a step-size h and calculate 


Y,, Y3,... from the second-order recurrence relation 


4hl, 54 hd, h 
Y,41 =(>-—_ | 7, ——_——]|Y,- ———— ](k,- 4k, +k 3 
r+1 (| ee (er | Sea + [| ee <3 os pad) ( ) 


where Y, = y(Xo), Y, has been previously obtained using some 
other numerical method and x, = Xo + rh. 


3. Y, is an approximation to y,. 


In the following example we use this procedure to illustrate the technique. 


Example 1 
Use Simpson’s method to solve the differential equation 


y=xt+y with y(0) = 0 (4) 
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using a step-size h = 0.1 for0 <x < 1. Initially generate the value of Y, at x = 0.1 
using a Taylor series method of appropriate order. 


Solution 
Since Simpson’s method has a local truncation error which involves h? it is 
worthwhile computing Y, very accurately. The Taylor series method of order n 
uses the recurrence relation 
h2 
YY, =i + eee 


Y" + ag yy 

where the derivatives may be obtained using the differential equation (4). We have 
Ti = 4, +7, Y,=1+Y', Y. =i, and so on. 

With h = 0.1 we have 


0.1)? 0.1)? 0.1)" 
Y,=¥, +01Y; ae a = Ye" +++ ¢ OF 


where Y =0, Yo=xo+ % =O, Yo=1+Y,=1, Y$.=Y$=1 andso on. 


O1y ay ty 0.1)" 
Hore ¥ 0404504 | Eas ee. ) 


rs” 


2 6 24 n! 
= 0.005 + 0.00016667 + 0.00000417 + 0.00000008 + --- 
= 0.00517092. 


The appropriate Taylor series method for eight decimal place accuracy for Y, is of 
order 5, 


With this value for Y, we use the procedure for Simpson’s method. We have 
lxj=4 
and Kix} = x. 
Hence the recurrence relation (3) becomes 
la = fea] Y,+ fea Y,-1 + [alte + 4x, + X,41). 
With h = 0.1, and noting that 
Npuy F Xe = oe = ee 


we have 


4 31 1 

39 z + 39 i + 579 (0-67). 

This second-order recurrence relation can be used to generate the solution using 
the initial conditions Y) = 0 and Y, = 0.00517092. In the following table the true 
solution 


Y,4+1 = 


y=e-x-1 


is also given for comparison. 


0 0 0 0 0 
1 0.00517092 0.0051709169 a1 % 10° 
Z 0.02140289 0.02140276 Soe 
3 0.04985897 0.04985881 1.6 x 107’ 
a 0.09182501 0.09182470 Eee | ame 
5 0.14872166 0.14872127 39 x 19% 
6 0.22211938 0.22211880 38x 10 
7 0.31375341 0.31375271 70x 10°? 
8 0.42554187 0.42554093 94x 10-’ 
’ 0.55960425 0.55960311 1.14 x 10°° 
10 0.71828328 0.71828183 1.42 x 10°° 


You may get slightly different 
results on your calculator. 
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All the results generated using Simpson’s method are correct to five decimal 
places. Comparing these results with the results obtained in Example 1 of Section 
1 for the same differential equation it can be seen that the results are much better 
for Simpson’s method than for Euler’s method, the trapezoidal method or the 
Taylor series method of order 2. Although quite a lot of work was necessary to 
compute a value for Y,, the results clearly make this effort worthwhile if great 
accuracy is required. 


Exercise 2 
Use Simpson’s method to solve the differential equation 


y’ = 3y + sinx with y(0) = 0. 


With h = 0.2, compute Y, and Y; given a value Y, = 0.0255. (Y,; was obtained using the 
trapezoidal method with a step-size h = 0.1.) 


Exercise 3 
Show that the use of Simpson’s method, with h = 0.1, applied to the differential equation 
y’ = 10y with y(0) = 1 


leads to the linear second-order recurrence relation 
Y4¢; @ 2 2 
By solving the auxiliary equation 
x? = 2x +2 
write down the general solution of the recurrence relation. 
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4.3 The numerical analysis of Simpson’s method 


As you saw in the last subsection, Simpson’s method does give very accurate 
results provided we obtain a good approximation to Y, using some other method. 
If hundreds or thousands of calculations are to be carried out then Simpson’s 
method is likely to give much better results with fewer calculations than any of the 
high-order Taylor series methods. 


However, there is one further drawback associated with Simpson’s method which 
we hinted at in Exercise 3: Simpson’s method cannot be used to solve certain 
problems. The use of Simpson’s method for a linear equation leads to a second- 
order recurrence relation. From Unit 1 we know that the general solution of a 
second-order recurrence relation contains two arbitrary constants. For example the 
general solution of the recurrence relation in Exercise 3 can be written as 


Y=Ax2 +8 xp 
where A = 2.7320508 and w = —0.73205081. 


In Unit 2 we saw that the general solution of a first-order recurrence relation 
involves only one arbitrary constant. For example the general solution for the 
differential equation 

y = iy 
10x 


iS y = Ae 


This general solution contains only one term while the general solution for 
Simpson’s method contains two terms. What is the correspondence between these 
two solutions? The general solution to the differential equation can be written as 


Vy aie Aer 
and 
Veni = Aes 
— 4e1r+h) 
ie Ae10*r x eioh 


pen 
=€ y;,- 
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Hence there is a recurrence relation for generating y,,, from y,. With h = 0.1 we 
have 

Yr+1 = Yr 

= 2.7182818y,. 

The general solution of this recurrence relation is 

Vn, = A(2.7182818)”. (5) 
Compare this with the solution of the recurrence relation 

Y, = A(2.7320508)" + B(—0.73205081)”. (6) 


There is a similarity between the first term in Y, and y,. The second term in Y,, 
B(—0.70205081)", is called a spurious solution because it is an extra solution 
introduced by the method. It does not correspond to anything in the original 
differential equation. How will this affect the results? 


We will apply Simpson’s method to the differential equation 
y = TOy with y(O) = 1. 


Suppose that we start the Simpson calculations by computing Y, = 2.718 using 
some other numerical method. Using this value and Y) = 1 we can calculate the 
appropriate values for A and B in (6). The particular solution is 


Y, = 0.9959 x (2.7320508)" + 0.0041 x (—0.73205081)". 


The particular solution for the differential equation is obtained by substituting the 
initial condition, yp = 1, into Equation (5) giving 


y, = (2.7182818)". 


Note that the first term in Y, is almost identical to the true solution, y,, while the 
spurious solution has a small coefficient (0.0041). As n increases the spurious 
solution gets smaller and smaller since (—0.73205081)”" tends to zero and hence its 
effect on the solution diminishes. We can conclude that in this case the spurious 
solution introduced by Simpson’s method does not seriously affect the results. 


In Exercise 4 you will see that, when Simpson’s method with h = 0.1 is applied to 
the differential equation 


y’ = -—10y with y(0) = 1, 


the spurious solution does seriously affect the results. Without going into the 
details of the analysis we state that Simpson’s method should not be used to solve 
differential equations of the form 


y’ = I(x)y + K(x) 


if I(x) is negative for any value of x in the interval under consideration. This is a 
generalization of the results obtained from Exercises 3 and 4. 


Exercise 4 
Use Simpson’s method to write down the recurrence relation for the differential equation 


y' = —10y 

with y(0) = 1. 

Using a step-size h = 0.1 show that the general solution of the recurrence relation is 
Y, = A(0.3660254)” + B(— 1.3660254)". 

What is the particular solution with Yo = 1 and Y, = 0.3679 ? 

Compare this with the true solution given by 
Yn = (0.36787944)”. 

What happens to Y, and y, as n increases? 


(You will have an opportunity to see the effect of this spurious solution when you use the 
computer package.) 


[Solution on p. 55] 
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Summary of Section 4 


1. The intervals of absolute stability for first-order methods are given in the 
following table. 


interval of 


method absolute stability 


Euler’s method 
trapezoidal method 

Taylor series method of order 2 
predictor-corrector method 


2. Simpson’s method can be used to generate approximate solutions to linear first- 
order differential equations using the procedure on page 40. 


3. Simpson’s method can give rise to very accurate solutions since the principal 
term in its local truncation error involves h°>. However, Simpson’s method 
should not be used to solve linear differential equations of the form 

y’ =I(x)y + k(x) 


if (x) is negative as, in this case, the spurious solution (the extra solution 
introduced by the method) will eventually swamp the solution. 
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5 The computer package NUMSOL 


5.1 Introduction 


This section is a description of the computer package NUMSOL. The package is 
designed to solve first-order differential equations using the methods described in 
this unit. It is intended that you should use this package at Summer School, 
although it will be available if you want to visit your local computer terminal. A 
general description of how to input data to packages is given in Unit 1 Section 6. 
This package is very similar to the package RECREL described in Unit 2. 


To help you to decide which method to use for a particular problem I have 
summarized the properties of each method in the following table: 


principal term 
in the 
method l.t.e involves advantages disadvantages 


simple to use; not very accurate; 
explicit method interval of absolute 
stability is (—2,0) 


Euler’s method 


Taylor series explicit method requires computation 
method of of y”. Interval of 
order 2 absolute stability 

is (—2, 0) 
predictor- explicit method interval of absolute 
corrector stability is (— 2,0) 
method (Euler- 
trapezoidal) 


method absolute stability 
is (— 00,0) 


explicit method; interval of absolute 
very accurate stability is 

(—2.51, 0)* Requires 
the computation of 


mt 


y. andy". 


Taylor series 
method of 
order 3 


very accurate implicit method; 
for some problems | should not be used 
for y’ = I(x)y + k(x) 
if I(x) is negative 


Simpson’s method 


trapezoidal interval of implicit method 


* not given in the text earlier. 


To illustrate how the package works we begin with an example showing how the 
logistic equation can be solved at the terminal. 


Example 1 
Solve the logistic equation 


y . = 
sana with y(O) = 100 


on the interval [0,2] using the predictor-corrector method with h = 0.1. 


y = 10y(1 
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Solution 


The following terminal dialogue shows how I input this problem, after logging in 
and obtaining the course library program NUMSOL. Each of my responses is 
underlined to indicate that this is what I typed. The other characters in the 
dialogue were output by the computer. 


Option 10 specifies the differential 


OPTION? 10 
foe equation to be solved 


"= 2?10*«Y*(1 — Y/1000 


OPTION? 23 This option specifies the method to be used 


PREDICTOR-CORRECTOR METHOD (EULER-TRAPEZOIDAL) 
OPTION? 30 
INITIAL CONDITION Y(X(0)) = ?100 


OPTION? 32 


These two options specify the 
range of x values 


INITIAL X VALUE X(0) = 70 
OPTION? 33 


FINAL X VALUE X(N) =?2 


OPTION? 34 


STEP-SIZE H = ?0.1 


OPTION? SOLVE Having input all the data, you can Solve the problem 


5.2 The options for NUMSOL 


Before using the package at a terminal you should plan which options you are 
going to use to solve your problems. A description of each option available for the 
package is given below. 


Command options 


OPTIONS — this tells the program to print a list of the available options. 


SOLVE -—this tells the program to run the problem after all the data has been 
input. If some data is missing an error message will be printed. 

HELP —this tells the program that you are stuck and need advice. 

LIST —this tells the program to print the problem and data previously input 


so that you can check the problem you are solving. 
STOP —this is the only way to stop the program. 


ANSWER —this gives the solution to the computing exercises so that you can 
check that your answer is correct. To the question 


EXERCISE? 


you respond with the number of the exercise you are interested in. 


Problem options 
10 —Input m(x, y) 
To the question 
Y=? 


, ae Valid input expressions are 
you respond with a valid input expression for m(x, y). described in Unit 1. 
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For linear equations you must input 
I(x) x y + k(x) 


in that order so that the package can recognize the linearity. 


11 —Order and derivatives for Taylor series method 
To the question 
ORDER = ? 


you respond with the order (<5) of the Taylor series method. When 
you have specified the order the program will ask you for each of 
the higher derivatives that the method requires. 


For example if you want to use a Taylor series method of order 2 
for the differential equation 


y = ayes 
you would respond to ORDER = ? with 2. The package would then 
ask you for the second derivative 

Y" =? 


to which you would respond Y + X*Y’ 


(The differential equation would have been input using Option 10 so 
Y’ is already known.) 


i —Change order of Taylor series method 


If you wish to change the order of the Taylor series method without 
having to input previous derivatives again, you would use this 
option. To the question 


ORDER = ? 


you would respond as in Option 11 and if further derivatives are 
required you will be asked for them. 


For example if you use the same differential equation as in the 
example in Option 11 but wish to change to the Taylor series 
method of order 3, to the question ORDER = ? you would respond 
with 3. The package already has Y’ (Option 10) and Y” (Option 11) 
and so only Y” is missing. The package will ask for this and to the 


question 
Note that only single quotes, 
YY =@ repeated as necessary, are 
used for derivatives e.g. Y’’ not 
you would respond2 * Y' + X*Y"’ ae 


Method options 

20 —Euler’s method 
This uses a first-order recurrence relation to solve the differential 
equation. The initial y value is input using Option 30. 

21 —Taylor series method 


The order of the Taylor series method and any further derivatives 
required are specified using Option 11 and the order can be changed 
using Option 12. The method uses a first-order recurrence relation to 
solve the differential equation. 


22 —Trapezoidal method 


This method can only be used to solve linear differential equations 
using a first-order recurrence relation. 
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23 —Predictor-corrector method 


This method uses Euler’s method to predict and the trapezoidal 
method to correct. Only one initial y value needs to be specified 
(Option 30). 

24 —Simpson’s method 


This method can only be used to solve linear differential equations. 
The method involves the use of a second-order recurrence relation 
for which two initial y values are required (Options 30 and 31). Y; 
may be computed using a one-step method with small step-size if it 
is not known. 

Data options 

30 —Jnitial condition Yo 


This specifies the initial y value. 


31 —Jnitial condition Y, 


This specifies the second initial condition for Simpson’s method. 


32 ——Xo_ 


—This is used to input the initial x value. 


33 —Xn 


This is used to input the final x value. 


34 —The step-size h 


This parameter specifies the step-size for the numerical method. 


Print options 
40 —QOutline printout 


If you only want selected output you should use this option to 
specify k such that every kth term in the sequence is printed. 


41 —Full printout (default) 


All terms in the sequence will be printed. 


42 —Print solution only 


Only the last term in the sequence is printed. (Intermediate terms 
will not be printed.) 


5.3 Computer exercises 


The following exercises will require the use of the computer package. Before you 
run the exercises at a terminal you should plan how you are going to tackle the 
problems you wish to solve. Judicious preliminary work will enable you to make 
the best use of your limited time at the terminal. You may not have time to solve 
all these problems. 


Exercise 1: This exercise will reproduce the results of the television programme. 
Use the package to solve the logistic equation 
v= 10 = TT with y(0) = 100 
forO <x <1, 
by each of the following methods 


(i) Euler’s method, 
(ii) predictor-corrector method, 
(iii) Taylor series method of order 2. 


A default option is the one the 
program will use if no other is 
specified. 


MST204 19.6 49 


For each method use a step-size h = 0.1. 


If you have time, look at what happens for larger step-sizes. 


Exercise 2: To investigate stability 


The differential equation 


1 
y ge ae with y(0) = 1 


is to be solved on the interval [0, 10] using each of the following methods 


(i) Euler’s method 
(ii) trapezoidal method 
(111) Simpson’s method 


with step-sizes h = 1, 0.5 and 0.1. 
(To reduce output, print the results only at the integer values of x.) 


Do the results agree with the stability analysis in Sections 3 and 4? 


Exercise 3: To investigate accuracy 
What methods can be used to solve the differential equation 
y =-2 with y(0) = 1 
on the interval [0,1]? 
Compare the numerical results obtained using Taylor series methods of order 2 and order 3 
and a step-size of h = 0.1. To get more accurate results use h = 0.05. 


d 
Hint: —e” = y’e” = —e?’, 
dx 


Is it more accurate, for this problem, to use a step-size h = 0.05 with a Taylor series method 
of order 2 or to use a step-size h = 0.1 with a Taylor series method or order 3? 


Exercise 4: To investigate stability 
Experiment with the differential equation 
y = —30y with y(0) = 1 


for 0 < x <1, to test the stability of the methods in the package for different values of h. 
The stability of Euler’s method for this equation is discussed in Section 3. 


Exercise 5: To investigate stability 

Use Simpson’s method with h = 0.1 to solve the differential equation 
y’ = —10y with y(0) = 1 

forO<x <2. 

Confirm the results of Exercise 4 in Section 4. 

(Take Y, = 0.3679 and then Y, = 0.3660.) 


[Solutions to the computer exercises are not given in the unit. A check solution can be 
obtained at the terminal using the option ANSWER. | 


6 End of unit exercises 


Exercise 1 


Compute the first 3 terms in the sequence if the trapezoidal method is to be used with a 
step-size h = 0.1 to solve the differential equation : 


y =y+2x with y(0) = 1. 

Exercise 2 

The predictor-corrector method is to be used with h = 0.1 to solve the equation 
y = ty with y(0) = 1. 


Write down a single recurrence relation which can be used to generate the sequence. What 
is the particular solution for this recurrence relation? 


The valid expression for e” is 


EXP(Y) 


50 


Exercise 3 
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Is the recurrence relation Y,,, = Y, + 5(4Y, + Y),,) consistent with the differential 


equation 


y’ = m(x, y) ? 


Would you expect this method to give better results than the trapezoidal method for small 


values of h? 


Exercise 4 


For what values of ha is the method in Exercise 3 absolutely stable when applied to the 


problem 
y= ay 
with y(0) = 1? 


Exercise 5 


Simpson’s method is to be used to solve the linear differential equation 


y’ = 3y + 2x 
with y(0) = 1. 


Derive the recurrence relation which can be used to generate the solution for a given step- 


size h. 


What other information do you require and how would you get it? 


[Solutions on pp. 55-56 | 


Appendix 


Solutions to the exercises in Section 1 
1. The trapezoidal method utilizes the recurrence 
relation 


h 
Yp+1 = q; ~ iad < I 44h 


For the differential equation 
y’ = 3y + sinx 


we have 
h 
Y,41= Y, + ces + sinx,) + (3Y,4, + sinx,4,)]. 
Collecting up the terms gives 
3h 3h a 
Y,41,1— a ie Y,|1 + > Se 5 (sin, + sinx,+1) 


i.e. 


(1 + 3h/2)Y, + (h/2)(sinx, + sinx,41) 
1 — 3h/2 


Ye+1 st 


[Alternatively we could just use the procedure box to write 
down the recurrence relation directly. | 


With h = 0.2 we have 


1.3Y, + 0.1(sinx, + sinx,4,) 
0.7 


Peat — 


The figures in the following table are quoted to 5 decimal 
places. 


0 | 0.02838 | 0.13672 | 0.39020 


0.02460 | 0.12308 | 0.35304 
po 0.00378 | 0.01364 | 0.03716 


Whilst the results are not particularly accurate, because of 
the rather large step-size, the ‘shape’ of the approximate 
solution is correct. The maximum error is approximately 0.04 
at x = 0.6. 


2. The formula for the trapezoidal method is 
h / / 
tpet — f 5 ate ss b er 1). 


For the non-linear differential equation 
y’ = siny 


we have 
ne : 
¥.44= ¥, + oe Y, +:sin Y,4 1). 
Collecting the terms in Y,,, and Y, we have 


ie, i . 
| rete =. + sett 
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Putting h = 0.1 gives 
Y,+1 — 0.05sinY,,, = Y, + 0.05sin Y,. 

With y(0) = Yo = 1, to get Y, we need to solve the equation 
Y, — 0.05 sin Y,; = 1.04207355. 


The solution cannot be written down immediately and the 
Newton-Raphson method, as described in Unit 18, Section 2, 
would have to be used to determine Y,. This method 
involves a considerable amount of extra work at each step. 


3. Assuming that Y, lies on a solution curve so that 


Y,=y,, Y,=y, and Y;,=y, 


we have, for the Taylor series method of order 2, 
2 


Y41=)r thy, + a Yr: 


The true solution y,;;, can be determined using a Taylor 
series expansion at x, as 


h? : ae 
fet = ee te te ee 
Hence the local truncation error is 
3 4 
14 Pe SS eur aie Sees 


The principal term in the local truncation error is 


3 


a es 


4. Assuming that Y, lies on a solution curve so that 


yr YY," Y,=¥, a Fey, 
we have, for the Taylor series method of order 3, 
h2 3 
YT, = r + h : + aay : + as a 
+41 =) y 7 y 6 y 


The true solution satisfies 


: 2 3 h* ho 
+1 = Vee: eee Pee ot 


+: 


Hence the local truncation error is 


Y,41 —-) +1 = 2 oo 
: : 34° ae 


The principal term in the local truncation error is 
4 
= ho, 
24° 
If the step-size is halved so that h* = h/2 the principal term 
in the local truncation error becomes 
(h*)* (h/2)* h* 


(iv) _ p62 


24°" 14°" 384 


The local truncation error has been reduced by a factor of 
approximately 16. Note however that about twice as many 
calculations would be needed to obtain an approximation at 
a particular value of x so that we would not expect the 
accumulated error at this value of x to be reduced by a 
factor of 16. 


5. The Taylor series method of order n is given by 


h2 
Y41=Y,+h¥,+—¥74+-: 


h" 
‘. ees y”, 
2 e n! 
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Assuming that Y, lies on a solution curve so that 


Ty FS ek TP ee 
we have 
h? h" 
Fee toe ee 
Z n! 
The true solution satisfies 
h2 
Vr+i at 57 t os 
h” n+1 yr*2 
Ee oe (n) - epeeenineene aid 1) (22) 
a ai Pe a 
4 WBS 
Hence the local truncation error is 
yn*i n+2 
ee er (n+ ee (n+ 2) 
pe ee (n+2)!° 


6. Here is the table of errors, for h = 0.1, using the Taylor 
series method of order 2 to solve the differential equation 


mx +y 
with y(0) = 1. 


0.005171 
0.021403 
0.049859 


— 0.000171 
— 0.000378 
— 0.000627 


0.005 
0.021025 
0.049232 


— 0.000044 
— 0.000098 
— 0.000163 


0.005171 
0.021403 
0.049859 


0.005127 
0.021305 
0.049696 


At x = 0.1 the errors are —0.000171 and —0.000044. Hence 
the error has been reduced by a factor of 171/44 ~ 3.89. 
Similarly at x = 0.2 and 0.3 the errors have been reduced by 
factors of approximately 3.86 and 3.85 respectively. 


We can conclude that, for the Taylor series method of order 
2, halving the step-size reduces the errors at particular values 
of x by a factor of approximately 4. (Note that if we halve 
the step-size the local truncation error is reduced by a factor 
of approximately 8.) 


7. Assuming that the values on the right-hand side are 
correct, we have 


Yr41=)r + hyyst. 
The Taylor series expansion for y,,, is given by 
h? ee 
te = Pl eC 


h3 
: Yra1 = Vr + hy, + heyy + yr + Ss 


The Taylor series expansion for y,,, is 


2 3 


 f h ” h wm 
Vr+1 =r tie ed ol ae 
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Hence 
2 3 
7; Suier = y, + yy, FO 
+1 — Yr+1 7 y 3 y 
The principal term in the local truncation error is 
h2 
ar ° 


This term involves h*, as does Euler’s method. 


8. (i) The principal term in the local truncation error for 
the trapezoidal method is 


3 


ID": 
If we halve the step-size, putting h* = h/2, the principal term 
in the local truncation error becomes 

as Boe 

ron se x et: 
Hence the local truncation error will reduce by a factor of 


approximately 8. 


The Taylor series method of order 2 also has a local 
truncation error whose principal term involves h? and so 
halving the step-size would reduce the local truncation error 
by a factor of approximately 8 for this method too. 


(ii) From the results of Exercise 6, for the Taylor series 
method of order 2, we would expect the accumulated errors 
for the trapezoidal method to reduce by a factor of 
approximately 4 if the step-size were halved. 


9. At x=0.1 the error has been reduced by a factor of 


92/23 = 4.00. Similarly at x = 0.2 and x = 0.3 the errors have 


been reduced by factors of 4.00 and 4.02 respectively. 


We conclude that, for this problem, halving the step-size 
reduces the errors at particular values of x by a factor of 
approximately 4 for the trapezoidal method. 


Solutions to the exercises in Section 2 


1. Separating the variables (see Unit 2) we have 


/ 


y S 
yal — y/l000) = ©) 


To get the left-hand side into the required form we use 
partial fractions to write 


1 .. . B 
y(1—y/1000) y (1 — y/1000) 
This leads to the two equations for A and B as 


tak 


se 
ee 
Hence Equation (1) becomes 


y’ y’ y’ 
a ee ee ee 
y(i —y/1000) y 1000—y 


Integrating both sides of Equation (1) gives 
logy — log(1000 — y) = 10x + C, 
1.€. 


y 


22 De iy eC 
ie] = 


oe| 
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giving 
y 
—_-__ = Ae!™, 
1000 — y 
This simplifies to 


_ 1000 Ae*™ 
ee 


For the initial condition y(0) = 100 we have 


ba 


1 
Hence A = -. 
9 


The particular solution is 
1000e'°* 
et ee 


1000 
On 1ST 


i 


Le. y= 


Here is a sketch of this solution. 


Sette 
pagdeanecezes 
ee 


ttt 
seagenens 


2. Using the predictor-corrector method to evaluate Y; and 
Y, we proceed through the four steps of the algorithm: 


1000 
(using the differential equation). 


(ii) Predict Y* = Y, + hY5 = 672.52255 


Y. 
(i) Evaluate Y5 = 10¥,(1 ee | = 2447.79 


(Euler’s method). 


= 2202.3597 


= Y% 
(iii) Evaluate Y$' = 10Y3}1 — : 
1000 


(using the differential equation). 


h 
(iv) Correct Y3 = Y2 + ae + Y¥') = 660.25103 
(trapezoidal method). 


Hence our approximation at x=0.3 is Y3 = 660.25103. 
Repeating the steps for Y, we have: 


Y3 
Evaluate Y3; = 10Y3{1 — 
(i) Evaluate Y3 | 1000 


| = 2243.1961 


(ii) Predict Yt = Y, + hY4 = 884.57064. 


* 


YZ 
= 1021.0542. 
ia 1.0542 


(iii) Evaluate YZ’ = 10¥e(1 — 


h 
(iv) Correct Y, = Y3 + ats + Y#') = 823.46355. 


Hence our approximation at x = 0.4 is Y, = 823.46355. 
3. For the differential equation 
with y(0) = 1 


and h = 0.1, the predictor-corrector method gives: 


y' =logy+x 
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(i) Evaluate Yo = log Yo + xo = logl + 0 =0. 
(u) Predict Y¥ = Y, + hYo = 1. 
(11) Evaluate Y*' = log Y* + x, = 0.1. 


h 
(iv) Correct Y, = Yo + 5! Yo + Y?’) = 1.005. 


Hence Y; == 1.005 at X4 — 0.1. 


(i) Evaluate Y, = log Y, + x, = 0.10498754. 
(11) Predict Y$ = Y, + hY{, = 1.0154988. 
(iti) Evaluate Y#' = log Y¥ + x = 0.21537988. 


h 
(iv) Correct Y, = Y, + 5% + Y3’) = 1.0210184. 


Hence Y, = 1.021 at x = 0.2. 
4. For the differential equation 
y =siny 


with y(0) = 1 and h = 0.2 the predictor-corrector method 
gives the following approximations at x = 0.2 and x = 0.4; 


(i) Evaluate Yo = sin Y, = 0.84147099. 
(ul) Predict Y¥ = Yo + hYo = 1.1682942. 
(11) Evaluate Y*' = sin Y* = 0.92008374. 


h 
(iv) Correct Y; = Yo + 3 (Yo + Y#’) = 1.1761555. 


Hence ¥Y, = Livéata = 02. 


(i) Evaluate Y;, = sinY, = 0.92313471. 
(i) Predict Y¥ = Y, + hY, = 1.3607824. 
(11) Evaluate Y%' = sin Y¥ = 0.97802802. 


h 
(iv) Correct Y, = Y; + ioe + Y3') = 1.3662717. 


Hence Y, = 1.366 at x = 0.4. 


3. Y*¥s1= Y,+hY,= Y, + hm(x,, Y,) using the differential 


equation. 


h 
Yat = ¥, 5s at = ye 


h 
= Lp 5g > [m(x,, Y,) + M(Xp415 oa i)]. 
Substituting for Y*,, in this second equation gives 
h 
Ae S [m(x,, Y,) + m(x,, Y, + hm(x,, Y,))] 


As this does not have Y;,, on the right-hand side the 
predictor-corrector method is explicit. 


Solutions to the exercises in Section 3 
1 
I. (i) P(x,, z,, | oe 1,h) = 3h + ea ores 1)- 
Hence 
1 , / 
P(X+. Vrs Vr, 9) — 3 a 2y;) 


= y, = M(x,,y,). 


Thus the method is consistent with the differential equation. 


1 
(ii) P(X, tn Y,+1,h) _ GY, se gre: 
Hence 


1 
P(X +, Vrs Vr, 0) si grr i yy) 


1 1 
= 5 = yr, Vr). 
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Thus the method is not consistent with the differential 
equation. 


2. (X,, Y,, Y-41,h) = aY, + bY, 41. 
Hence 
P(Xr5Vr,Vr,0) = ay, + by; 
= (a + b)m(x,,y,). 
The method is consistent with the differential equation 
y' = m(x,y) 
only ifa+b=1. 


If a = b = 4 we have the trapezoidal method. If a = 1, b= 0 
we have Euler’s method. 


3. From Exercise 2 with a = b = 4 we know that the 
trapezoidal method is consistent. From Theorem 1 the 
method is also convergent. Since the principal term in the 
local truncation error for the trapezoidal method involves h?, 
Theorem | tells us that the global error at x* is of the form 


Yy + — y(x*) ~ Ch? 
for small values of h. Thus halving the step-size would reduce 
the global error at x* by a factor of approximately 4. 
4. From Exercise 5 in Section 2 we have 


h 
Y.j,=¥4 = [m(x,, ¥,) + m(x,, Y, + hm(x,, Y,))]. 


Hence 
P(X, Ys Ye+1, h) = z[m(x,, ¥,) + m(x,, ¥, + hm(x,, Y,))]. 
Thus 
P(Xr,Vr, Vr, 0) = 5 [m(x,,y,) + m(x,,y,)] 
= m(x,, yr). 


The predictor-corrector method is consistent with the 
differential equation. 


5. Euler’s method uses the recurrence relation 
Y-41= Y, +hyY,. 
For the differential equation 
, = oy 
with y(0) = 1 and h = 0.1 we have 
Y,+1 = Y, — 30 x 0.1Y, 
s —2Y,. 


With Y) = 1 we have the particular solution for this 
recurrence relation as 


Y, = (—2)" 


This is an increasing oscillating solution. (The sequence 
generated is 1, —2, 4, —8, 16, —32, 64...) 


The true solution is y = e~ *°* which tends rapidly to zero. 
Something has clearly gone wrong with the method. 


6. Euler’s method is absolutely stable if —2 < ha <0. 
Hence for the differential equation 


y= 10 
a = —100 giving 
—2< —100h <0 


i.e. 100h < 2. (The inequality —100h < 0 is always satisfied 
for positive h.) 


1.e. h < 0.02. 


7. The differential equation 


/ 


y'=—xy+ 


ee 


is linear with I(x) = —x and k(x) = vee 
x 


The condition for absolute stability in Euler’s method is that 
—2<I(x)h <0. 

The minimum value for I(x) is —10 since 0 <x < 10. 

Hence we require 
—2 < —10h. 

Le. Wh = 2 


h < 0.2. 


8. The theory tells us that y lies between 1000 and 2000 in 
the solution of the logistic equation 


y : 
ae — —— h y(0) = 2000. 
y 101 eal with y(O) 


Hence the minimum value of 


d 
ar eee 
dy 50 
3. ae 
occurs when y = 2000 giving o =: — 30. 
We require 
—2< qot < 0. 
dy 
Le... Ae 
h< S 
15 


Thus, provided h < 1/15, Euler’s method is absolutely stable. 
9. This is a differential equation of the form 

y' =m(x,y). 
At x, we write 


m,(y) = m(x,, Yr). 


6x 
In this case m(x, y) = — 
y 


so that 
6x, 
m,(y) = —. 
y 

Now 
dm, _—‘6X, 
7 


The condition for absolute stability in Euler’s method is 


dm, 


—2<h < 0. 
dy 


x lies between 0 and 2 whilst y lies between 1 and 5. Hence 


= dm, 
the minimum value of occurs when x, = 2 and y= 1 
y 


giving 
= > —12. 

Hence the method is certainly absolutely stable if 
—2 < —12h. 
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1 
Le Rs, 

6 
(In reality we could get a much less restrictive bound on h by 
obtaining the solution using the separation of variables 
method and then finding the minimum value of — 6x/y’, but 
that would defeat the object of the exercise!) 


Solutions to the exercises in Section 4 


1. For the Taylor series method of order 2 we have 


h? 
Y,41= Y, +hY, + Yr. (1) 


Applied to the problem 
y =a 


we have Y= aY, and Y, = aY) = «’Y, (differentiating the 
differential equation). Substituting into Equation (1) above 
gives 


pe | 
ad = Y, + haY, + 


Y, 


h*a? 
= (14+ ha + 5 )y. 


This is a recurrence relation of the form 
a Sg ay, 


where 
2.2 


h*a 
a=1+ha+ ee 


Note that this is the same equation for a that we had in the 
predictor-corrector method in tape frame 3. Hence the graph 
will be the same as in tape frame 4 and consequently the 
interval of absolute stability (|a| < 1) is (—2,0). 


In the logistic equation, with y(0) = 500 the solution 
d 
increases rapidly to 1000 and 2 lies between 0 and — 10 
- 


where 


a 
1000) 
As the interval of absolute stability for this method is (—2,0) 


m(y) = 101 — 


d 
the method is stable provided h ie Sn hche he ee 
y 


2. Using the procedure for Simpson’s method applied to 
linear differential equations with I(x) = 3, k(x) = sinx and 
h = 0.2 we obtain the recurrence relation 


2.4 3.6 


0.2 
a aat Y,-, += (sinx,_, + 4sinx, + sinx,+1) 


33 2.4 
ae 
= 3, + hodjes © eee + 4sinx, + sin x,+1) 


Using the initial values Y) = 0 at x» =0 and Y, = 0.0255 at 
x, = 0.2 we have: 


1 
Y, = 00255 + 15 x 0+ oT + 4sin0.2 + sin0.4) 


= 0.12417464; 
Y, = 0.12417464 + 1.5 x 0.0255 


1 
+ D (sin 0.2 + 4sin 0.4 + sin 0.6) 


= 0.35584007. 
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3. Using the procedure for linear equations with I(x) = 
k(x) = 0 and h = 0.1 we obtain the recurrence relation as 


4x 0.1 x 10 3+ 0.1 x 10 
3 —0.1 x 10 3 — 0.1 x 10 


— re & + Fy Somer ee 


| }¥,1 +0 


This is the required recurrence relation. 
The auxiliary equation (cf. Unit 1, Section 2) is 
MO SON? of. x? = 2x — 2 SU. 


The roots of this quadratic are 


so 4 op / 5. 
Thus 4 = 1 + ./3 = 2.7320508 and 
uw =1— ./3 = —0.73205081. 
The general solution for this recurrence relation is 
Y,=AxA"+Bx p" 
where 4 = 2.7320508 and pw = —0.73205081. 


4. Again using the procedure for linear equations with 
I(x) = —10, k(x) = 0 and h = 0.1 we have 
4x 0.1 x (—10) y 
3 — 0.1 x (—10) 
3 + 0.1 x (—10) 
3 — 0.1 x (—10) 


as = ¥o+ $7,.1. 


ek = | 


}¥1 +0 


This is a linear, constant-coefficient, homogeneous, second- 
order recurrence relation. Its auxiliary equation is 


2 


x*=—-x+$ or x?+x-43=0. 


The roots of this equation are 


cast ae SPT ie oe 


A, b=; 
. 2 


-§ = 
oS = ~ 03660254 and 


Ss — — 1,3660254. 


Hence the general solution is 
Y, = A(0.3660254)" + B(— 1.3660254)”. 


If Yo = 1, and Y, = 0.3679 the particular solution is obtained 
from solving for A and B in 


Y=A+B=1 

Y, = 0.3660254A — 1.3660254B = 0.3679. 
This gives A = 1.0010823 and B = —0.0010823. 
Hence the particular solution is 


Y, = 1.0010823 x (0.3660254)” 
—0.0010823 x (—1.3660254)". 


Compare this with the true solution 
= (0.36787944)". 


The first term in Y, is almost exactly the same as the true 
solution while the spurious solution has a small coefficient 
(—0.0010823). However see what happens as n increases: 
(0.3660254)”" tends to zero while (— 1.3660254)” is oscillating 
and increasing. Inevitably the spurious solution will 
eventually dominate the values obtained for Y,. Below are 
the values I obtained using the recurrence relation with 

¥, = 1 and Y, = 0.3679. 


= 


Note that the solution decreases initially, but then the 
spurious solution takes over so that after x = 0.5 the solution 
is oscillating and increasing. 


Solutions to the end of unit exercises in 
Section 6 


1. The recurrence relation for the trapezoidal method is 


h 
ee = b a Str ~ Y 4 1) 


h 
= ie + —(Y, - 7.3 1) 


h 
2 blots tke a) 


Hence 
h h 
Yres(1 — 4 = (1 + 4 +(x, + X, + i): 
2 2 
(Alternatively we could have used the procedure for linear 
equations with /(x) = 1 and k(x) = 2x.) 
With h = 0.1 we have 


- LOSY, + 01x; + X,+1) 
0.95 


2. The predictor-corrector method applied to the 
differential equation 

yo Ty with y(0) = 1 andh=0.1 
gives 


(i) Evaluate Y, =7Y,. 
(ii) Predict Y* r+ = = PRY, = ¥, + hy, = 1.7Y,. 
(i111) Evaluate Y*)., = 7Y*, = = 11.9Y,. 


(iv) Correct Y,4, = Y, + Hc + Y*' ,). 
: 
= Y,+ 5(7Y, + 11.9Y,) = 1.945Y,. 


Hence Y,,, = 1.945Y,. 


The particular solution for this recurrence relation with 
Yo = | is 


Y, = (1.945)". 
3. For the recurrence relation 


h 
foe gun a; 2 54Y, + ris 1) 


we have 


P(X>, = ¥,+1,h) a 3(4Y;, - ee 1). 


Thus 
P(X,5 Vrs Vr50) = $(4y; + ye) 


oad yy = m(X,, Vr). 


Hence the method is consistent with the differential equation. 


To compare this method with the trapezoidal method we 
compute the local truncation error. As it is an implicit 
method we can only compute the principal term, by 
assuming that all right-hand side values are correct; 


: h pe ! 
1.€. T+ me 2 = 5 yr + Yr+1)- 
The Taylor series expansion for y,+, is 


h2 
OE ee Saat SS ore Pee 


2 
Hence 
: h? 
Y, =y), hy’, —y, +— ee 
41 =)rt ny, + 5 y 10” 
The Taylor series expansion for y,+; gives 
h2 h3 
r+1 = , + hy + —y, ty ot 
Yr+i=y - D y 6 - 
: Y a= 3h? ”" h° m 
ey (ag | Vr+i= 10°" 15)" 


2 


3h ae 
The principal term, — T97" involves h? while the principal 


term in the local truncation error for the trapezoidal method 
involves h?. We would thus expect better results using the 
trapezoidal method. 


4. For the problem 
y = ay 


we have 


h 
Vee = = + ee + at, 1). 
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Collecting up terms gives 
y ' ha\ Wiad 4ha\ 
r+1 5 pas: 5 ’ 
Le. 
(1 + 4ha/5) 
a Sa ees Se 
(1 — ha/5) 
For absolute stability we require 


! + Aha/5 


1. 
i 


and the method is 


: ae Aha ha 
If ha is positive |1 + =< so = 


absolutely unstable. 


4ha ha 10. 
If ha is negative |1 + os <jl— a provided ha > a in 


which case the method is absolutely stable. 
10 
Hence the method is absolutely stable if — = < ha < 0. 


5. Using the procedure for linear equations with I(x) = 3 
and k(x) = 2x we obtain the recurrence relation for 
Simpson’s method applied to the linear equation as 


oe 4h - pee 
ra al Re = 5 r fk P= 
2h 


= 3 30 -w + 4x, + Saat 


Since x,_1 + X,;+1 = 2x, this can be simplified to 


e 4h 1+h y 3 4h 
= |———_ ——] Y,- ———]x,. 
eee ee eee 
To get started with this second-order recurrence relation, we 
require a value for Y,. This can be obtained, for example, by 


using a one-step method with very small step-size to compute 
an accurate approximation to Y;. 
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